mRNA-1273 and Ad26.COV2.S vaccines protect against the B.1.621 variant of SARS-CoV-2

Background Since the emergence of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in 2019, viral variants with greater transmissibility or immune-evasion properties have arisen, which could jeopardize recently deployed vaccine- and antibody-based countermeasures. Methods Here, we evaluated in mice and hamsters the efficacy of a pre-clinical version of the Moderna mRNA vaccine (mRNA-1273) and the Johnson & Johnson recombinant adenoviral-vectored vaccine (Ad26.COV2.S) against the B.1.621 (Mu) variant of SARS-CoV-2, which contains spike mutations T95I, Y144S, Y145N, R346K, E484K, N501Y, D614G, P681H, and D950N. Findings Immunization of 129S2 and K18-human ACE2 transgenic mice with the mRNA-1273 vaccine protected against weight loss, lung infection, and lung pathology after challenge with the B.1.621 or WA1/2020 N501Y/D614G SARS-CoV-2 strain. Similarly, immunization of 129S2 mice and Syrian hamsters with a high dose of Ad26.COV2.S reduced lung infection after B.1.621 virus challenge. Conclusions Thus, immunity induced by the mRNA-1273 or Ad26.COV2.S vaccine can protect against the B.1.621 variant of SARS-CoV-2 in multiple animal models. Funding This study was supported by the NIH (R01 AI157155 and U01 AI151810), NIAID Centers of Excellence for Influenza Research and Response [CEIRR] contracts 75N93021C00014 and 75N93021C00016, and the Collaborative Influenza Vaccine Innovation Centers [CIVIC] contract 75N93019C00051. It was also supported, in part, by the National Institutes of Allergy and Infectious Diseases Center for Research on Influenza Pathogenesis (HHSN272201400008C) and the Japan Program for Infectious Diseases Research and Infrastructure (JP21wm0125002) from the Japan Agency for Medical Research and Development (AMED).

INTRODUCTION
Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the etiological agent of coronavirus disease 2019 (COVID- 19), has caused hundreds of millions of infections worldwide, with more than 5 million deaths. Vaccines targeting the SARS-CoV-2 spike protein were developed within 1 year of the start of the pandemic. Several of these (mRNA and adenoviral vectored) are remarkably effective in protecting against severe COVID-19, with efficacy rates ranging from 75%

Context and significance
New variants of SARS-CoV-2 have emerged that are more resistant to neutralizing antibodies. Prior to Omicron, the Mu variant (B.1.621) was one of the more neutralizingantibody-resistant variants of SARS-CoV-2. Here we evaluated the efficacy of two vaccines for protection against the Mu variant in mice and hamsters. Immunization with the mRNA-1273 or Ad26.COV2.S vaccine induced serum neutralizing antibody responses and protected the animals against disease after challenge with B.1.621. Thus, immunity induced by these vaccines can protect against the B.1.621 variant of SARS-CoV-2 in multiple animal models.
Med 3, 309-324, May 13, 2022 ª 2022 Elsevier Inc. 309 ll to 95% depending on the vaccine and age of the individual. [1][2][3] Vaccines also protect against infection, and likely transmission, albeit at lower rates of 50%-70%. 4,5 The emergence of several SARS-CoV-2 variants with amino-acid substitutions in the spike protein has jeopardized the efficacy of current vaccines to protect against infection and disease. These variants can be more transmissible and also evade serum neutralizing antibodies. As an example, while the Delta variant (B.1.617. 2) had no appreciable effect on vaccine efficacy against hospitalization, vaccine-mediated protection against infection was markedly reduced. [1][2][3][4][5] Thus, evaluating vaccine efficacy against emerging variants of SARS-CoV-2 is important for deciding when to administer booster shots and determining if and when mono-or multivalent vaccines with variant spike antigens are needed.

RESULTS
Immunogenicity and protection by Ad26.COV2.S vaccine against B.1.621 challenge in 129S2 mice Groups of 5-to 6-week-old male 129S2 mice were immunized once with 10 8 , 10 9 , or 10 10 virus particles of fresh or freeze-thawed Ad26.COV2.S vaccine ( Figure 1A). Serum was collected 21 and 115 days later, and antibody responses were evaluated by ELISA for spike (S)-specific antibody responses. As expected, serum from control mice that received PBS did not bind to the S protein by ELISA ( Figure 1B). In comparison, serum collected from mice 115 days after immunization with a 10 8 , 10 9 , or 10 10 dose of Ad26.COV2.S contained anti-S-specific antibodies with geometric mean titers (GMTs) of 1:21,087, 1:19,870, and 1:80,599, respectively ( Figure 1B). No differences in GMTs were observed between animals immunized with the fresh or freezethawed Ad26.COV2.S vaccine (Figures S1A and S1B). The serum anti-S antibody response was higher in animals that received the 10 10 dose (p < 0.001; Figure 1B) than those immunized with the 10 8 or 10 9 dose of vaccine. A comparison of the anti-S response between 21 and 115 days after immunization revealed an approximately 3-fold increase in anti-S response over time ( Figure S1C). Serum samples also were tested for neutralization of SARS-CoV-2 by focus reduction neutralization test (FRNT). Whereas serum from the control animals did not neutralize WA1/2020 (B) Serum anti-S protein antibody response (EC 50 ) in control mice (black symbols) and mice immunized with 10 8 (red symbols), 10 9 (purple symbols), or 10 10 (blue symbols) of Ad26.COV2.S (****p < 0.0001, ***p < 0.001 by non-parametric one-way ANOVA with a Dunn's post-test). (C) Serum neutralizing titer (IC 50 ) against WA1/2020 N501Y/D614G (circles) or B.1.621 (triangles) measured by FRNT from 129S2 mice immunized once with 10 8 (red symbols), 10 9 (purple symbols), or 10 10 (blue symbols) of Ad26.COV2.S (****p < 0.0001, ***p < 0.001, ns = not significant by non-parametric one-way ANOVA with a Dunn's post-test). (D-F) 129S2 mice were challenged with 10 5 PFU of the WA1/2020 N501Y/D614G (circles) or B.1.621 (triangles) variant of SARS-CoV-2, and nasal washes (D) and lungs (E and F) were collected for analysis of viral RNA levels by qRT-PCR (D) and infectious virus by plaque assay (E and F) (****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05, ns = not significant by one-way ANOVA with a Dunnett's [D and E] or Dunn's [F] post-test). (B-F) Bars indicate the geometric mean values, and dotted lines are the limits of detection (LODs) of the assays. Animals at the limit of detection are arbitrarily assigned this value. These values are combined with those having values above the limit to determine the GMT. The results are from two independent experiments, and each symbol represents an individual animal. See also Figure S1.
Immunogenicity and protection by mRNA-1273 against B.1.621 challenge in K18-hACE2 mice Next, we evaluated the efficacy of pre-clinical non-good manufacturing practice (GMP) lots of the Moderna mRNA-1273 vaccine encoding a sequence-optimized 2 proline-stabilized S protein of Wuhan-1 for protection against B.1.621 in K18-hACE2 transgenic mice. Groups of 7-to 8-week-old female mice were immunized and boosted via an intramuscular route with a low (0.25 mg) or high (5 mg) dose of the mRNA-1273 or a control mRNA vaccine (mRNA-control; Figure 2A). Serum was collected 21 days after the second immunization, and inhibitory antibody responses were evaluated by FRNT against WA1/2020 N501Y/D614G and B.1.621. As expected, serum from mice immunized with the control vaccine did not inhibit virus infection ( Figure 2B). In contrast, serum from mice immunized with a 0.25 mg dose of mRNA-1273 neutralized the infectious virus with GMTs of 1:1,125 and 1:434 for WA1/2020 N501Y/D614G and B.1.621 viruses, respectively. Serum from mice immunized with 5 mg mRNA-1273 showed greater neutralizing activity against WA1/2020 N501Y/D614G and B.1.621, with GMTs of 1:19,751 and 1:15,130, respectively. A statistical difference in serum GMTs was observed between WA1/ 2020 N501Y/D614G and B.1.621 for the low dose (p < 0.05), but not the high dose (p > 0.5), vaccine. However, a pairwise comparison showed a 1.3-to 2.3-fold difference in serum neutralizing antibody titers between WA1/2020 N501Y/ D614G and B.1.621 for the high-(p < 0.05) and low-(p < 0.01) dose mRNA-1273immunized mice ( Figure S1E). We next evaluated the protective effect of the mRNA-1273 vaccine in K18-hACE2 mice. Mice immunized with mRNA-1273 or control mRNA vaccine were challenged via an intranasal route with 10 3 PFU of WA1/2020 N501Y/D614G or B.1.621 virus, and body weight was recorded for 7 days before a nasal wash and the left lung lobe were collected for viral-burden analysis. In the control groups challenged with the WA1/2020 N501Y/D614G virus, substantial weight loss was observed at 6 and 7 dpi. Immunization with 0.25 or 5 mg of mRNA-1273 significantly prevented weight loss (p < 0.05 and 0.01, respectively; Figure 2C). Similar protection against weight loss after challenge with B.1.621 was observed in the groups immunized with 0.25 (p = 0.19) or 5 mg (p < 0.05) of mRNA-1273 ( Figure 2D).
We quantified the amount of viral RNA in the nasal wash at 7 dpi. In control-vaccinated animals challenged with WA1/2020 N501Y/D614G, we detected $10 5 copies of the SARS-CoV-2 N gene transcript/mL in the nasal wash ( Figure 2E). Immunization with mRNA-1273 reduced the levels of WA1/2020 N501Y/D614G RNA by 20-to 45-fold for the 0.25 (p < 0.05) and 5 mg (p < 0.01) doses, respectively ( Figure 2E). In the groups vaccinated with the control mRNA and challenged with the B.1.621 virus, we detected 10 5 -10 6 copies of the N gene transcript/mL in the nasal wash ( Figure 2E). Immunization with 0.25 (p = 0.1) or 5 mg (p < 0.001) of mRNA-1273 reduced the N gene copy number by $25and 200-fold, respectively ( Figure 2E).
The amount of virus in lung homogenates was also quantified. In the control group challenged with WA1/2020 N501Y/D614G, $10 6 copies of the N transcript/mg and 10 5 -10 6 PFU/g tissue were detected ( Figures 2F and 2G). Immunization with 0.25 mg of mRNA-1273 resulted in markedly reduced ($20,000-fold, p < 0.0001) levels of viral RNA and no detectable infectious virus titer in the lung ( Figures 2F and 2G). Similar results were seen in mice immunized with 5 mg of mRNA-1273 ( Figures 2F  and 2G). A second group of immunized animals was challenged with B.1.621. In the control-mRNA-vaccinated groups, we detected high levels of viral RNA (N gene, >10 6 copies/mg) and infectious virus (6 3 10 5 PFU/g) in lung homogenates  K18-hACE2 mice vaccinated with 0.25 or 5 mg of mRNA-1273 also had markedly reduced, if not absent, lung pathology at 7 dpi compared with the control mRNAvaccinated and challenged animals ( Figure S2). Lung sections obtained from mice immunized with the control mRNA vaccine and challenged with B.1.621 showed severe pneumonia characterized by immune cell infiltration, alveolar space consolidation, vascular congestion, and interstitial edema. In comparison, mice immunized with either dose of mRNA-1273 did not show evidence of lung pathology. Overall, these data indicate that the mRNA-1273 vaccine protects against the B.1.621 variant of SARS-CoV-2 in K18-hACE2 mice.
Immunogenicity and protection by mRNA-1273 vaccine against B.1.621 challenge in 129S2 mice To corroborate our findings, we also tested the mRNA-1273 vaccine in immunocompetent 129S2 mice. Groups of 7-to 8-week-old female mice were immunized and boosted via an intramuscular route with 0.25 or 5 mg of mRNA-1273 or control mRNA vaccine ( Figure 3A). Serum was collected 21 days after the second immunization, and antibody responses were evaluated by FRNT. As expected, serum from 129S2 mice immunized with the control mRNA vaccine did not neutralize virus infection ( Figure 3B). In contrast, serum from mice immunized with 5 mg of mRNA-1273 robustly neutralized infection with GMTs of 1:40,066 and 1:38,675 for WA1/2020 N501Y/D614G and B.1.621, respectively. Sera from mice immunized with the 0.25 mg dose of mRNA-1273 inhibited infection of WA1/2020 N501Y/D614G and B.1.621 to a lesser extent, with GMTs of 1:2,665 and 1:2,407, respectively. No significant difference (p > 0.05) in vaccine-induced GMTs were observed between the WA1/2020 N501Y/D614G and B.1.621 viruses ( Figure 3B). Also, no significant difference (p > 0.5) was observed in a pairwise comparison of the neutralization titers against WA1/2020 N501Y/D614G and B.1.621 in serum from 129S2 mice immunized with high-and low-dose mRNA-1273 vaccine ( Figure S1F). We next challenged immunized 129S2 mice with 10 5 PFU of WA1/2020 N501Y/ D614G or B.1.621 via an intranasal route. Weights were recorded for 4 days before the animals were necropsied, and nasal washes and the left lung lobe were collected for viral-burden analysis. Control animals challenged with the WA1/2020 N501Y/ D614G virus lost $10% of their starting body weight by 4 dpi ( Figure 3C). Immunization with 0.25 and 5 mg of mRNA-1273 reduced the weight loss (p < 0.001 and 0.01, respectively; Figure 3C). Inoculation of control-vaccinated 129S2 mice with B.1.621 (E-G) 129S2 mice were challenged with 10 5 PFU of WA1/2020 N501Y/D614G (red symbols) or B.1.621 (blue symbols), and nasal washes (E) and lungs (F and G) were evaluated for viral RNA levels by qRT-PCR (E) and infectious virus by plaque assay (F and G) (****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05, ns = not significant by Mann-Whitney test). (B and E-G) Bars indicate the geometric mean values, and dotted lines are the LODs of the assays. Animals at the limit of detection are arbitrarily assigned this value. These values are combined with those having values above the limit to determine the GMT. The results are from two independent experiments, and each symbol represents an individual animal. See also Figures S1 and S3. resulted in $7% weight loss at 4 dpi ( Figure 3D), and immunization with 0.25 or 5 mg of mRNA-1273 prevented this weight loss (p < 0.05).
In animals vaccinated with the control mRNA and challenged with the WA1/2020 N501Y/D614G virus, we detected $10 6 copies/mL of N gene transcript ( Figure 3E) in the nasal wash at 4 dpi. Immunization with 0.25 or 5 mg of mRNA-1273 reduced WA1/2020 N501Y/D614G viral RNA levels in the nasal wash ($200-fold, p < 0.001). In mice vaccinated with the control mRNA vaccine and challenged with B.1.621, we detected $10 6 N gene transcript copies/mL in the nasal wash (Figure 3E). Immunization with 0.25 or 5 mg of mRNA-1273 reduced the viral RNA levels substantially ($200-fold, p < 0.001; Figure 3E). Viral burden also was measured in lung homogenates of these same animals. In the control mRNA group challenged with the WA1/2020 N501Y/D614G virus, we detected $10 8 copies/mg of N gene transcript ( Figure 3F) and >10 6 PFU/g of infectious virus ( Figure 3G) in lung homogenates. Immunization with 0.25 mg of mRNA-1273 effectively reduced the viral RNA (10 6 -fold, p < 0.001; Figure 3F) and infectious virus levels (1,000-fold, p < 0.001; Figure 3G), although some animals showed breakthrough infection. Immunization with a 5 mg dose of mRNA-1273 also reduced the levels of viral RNA (10 6 -fold, p < 0.001) and infectious virus (1,000-fold, p < 0.001; Figures 3F and 3G). Importantly, no breakthrough infection was detected. A second group of immunized animals were challenged with the B.1.621 virus. In the control-mRNA-vaccinated groups, we detected in the lung homogenates $10 7 copies/mg of viral RNA transcript ( Figure 3F) and $10 6 PFU/g of infectious virus ( Figure 3E). Immunization with 0.25 or 5 mg doses of mRNA-1273 significantly reduced B.1.621 viral RNA (10 5 -fold, p < 0.0001) and infectious virus (1,000-fold, p < 0.0001) levels in the lung ( Figures 3F and 3G) to a similar degree after WA1/2020 N501Y/D614G challenge.
The reduction in viral burden in the lungs of mice immunized with 0.25 or 5 mg of mRNA-1273 and challenged with WA1/2020 N501Y/D614G or B.1.621 corresponded with an absence of lung pathology at 4 dpi ( Figure S3). Lung sections from mice immunized with a control mRNA and challenged with WA1/2020 N501Y/D614G or B.1.621 showed evidence of immune cell infiltration and mild tissue injury characterized by vascular congestion and interstitial edema. Immunization with 0.25 and 5 mg of mRNA-1273 prevented lung pathology after challenge with both the WA1/2020 N501Y/D614G and B.1.621 viruses. Overall, these data indicate that the mRNA-1273 vaccine protects against the B.1.621 variant of SARS-CoV-2 in non-transgenic immunocompetent 129S2 mice.
Immunogenicity and protection by Ad26.COV2.S vaccine against B.1.621 challenge in Syrian hamsters We next evaluated the efficacy of a single dose of the Ad26.COV2.S vaccine in Syrian hamsters. Syrian hamsters are naturally susceptible to SARS-CoV-2 and considered an excellent model for COVID-19. 19,21-24 Groups of 5-to 6-week-old male hamsters were immunized once with the 10 8 or 10 10 dose of the Ad26.COV2.S vaccine ( Figure 4A). Serum was collected 21 days later, and antibody responses were evaluated by ELISA and FRNT. As expected, serum from control hamsters that received PBS did not bind the viral S protein ( Figure 4B). However, sera collected from hamsters immunized with the 10 8 or 10 10 dose of Ad26.COV2.S contained anti-S antibodies with GMTs of 1:7,506, and 1:40,913, respectively ( Figure 4B). The serum anti-S antibody response was higher in animals receiving the 10 10 vaccine dose (p < 0.05; Figure 4B). Serum samples were tested for neutralization of SARS-CoV-2 by FRNT. Serum from hamsters immunized with the 10 8 or 10 10 Figure 4C). A pairwise comparison found a $3fold reduction in serum neutralizing antibody titers between WA1/2020 and B.1.621 for the high-(p < 0.0001) and low-(p < 0.001) dose Ad26-COV2.S-immunized hamsters ( Figure S1G).
Seventy days after immunization, the hamsters were challenged via an intranasal route with 10 3 PFU of WA1/2020 or B.1.621. Four days later, the animals were necropsied, and a nasal wash and the left lung lobe were collected for viral-burden analysis. In the nasal wash of control animals challenged with the WA1/2020 virus, we detected $10 6 copies of the N gene transcript/mL of nasal wash ( Figure 4D) and $10 4 PFU/mL of infectious virus ( Figure 4E). Immunization with the 10 8 or 10 10 dose of Ad26.COV2.S did not reduce the viral RNA levels in the nasal wash (Figure 4D), although the infectious virus titer was decreased by 16-and 46-fold (p > 0.05 and <0.05, respectively; Figure 4E). Upon challenge of hamsters with B.1.621, we detected $10 6 N gene copies/mL ( Figure 4D) and $10 3 PFU/mL of infectious virus ( Figure 4E) in the nasal wash. Immunization with the 10 8 or 10 10 dose of Ad26.COV.2 reduced the N gene copy number 19-fold (p > 0.05; Figure 4D). In comparison, the infectious virus titer was reduced by 2-and 17-fold (p > 0.05 and < 0.05 for the 10 8 and 10 10 doses, respectively; Figure 4E) We also measured the viral burden in lung homogenates of hamsters. In the control group challenged with WA1/2020, the N gene copy number was $10 7 copies/mg ( Figure 4F) and the infectious virus levels were $10 7 PFU/mL ( Figure 4G). Immunization with the 10 8 or 10 10 dose of Ad26.COV2.S reduced the viral RNA (11-and 360-fold, p > 0.05 and < 0.01) and infectious virus (66-and 4,800-fold, p = 0.06 and < 0.001) levels ( Figures 4F and 4G). After B.1.621 challenge of the control group, we detected $10 7 N gene copies/mg of tissue ( Figure 4F) and $10 7 PFU/mL of infectious virus in the lung ( Figure 4G). Immunization with the 10 8 or 10 10 dose of Ad26.COV2.S reduced the viral RNA (7-and 540-fold, p > 0.05 and < 0.05) and infectious virus (100-and 25,500-fold, p = 0.13 and < 0.001) levels ( Figures 4F and 4G). A comparison between neutralization titers (half maximal effective concentration [EC 50 ] values) and infectious virus titers in the lungs in animals challenged with WA1/2020 revealed a significant correlation (R 2 = 0.67, p = 0.0012; Figure S4). In contrast, this association was not clearly observed in hamsters challenged with B.1.621 (R 2 = 0.24, p = 0.11; Figure S4), suggesting that other immune factors (non-neutralizing antibodies, T cells, or anamnestic B cell responses) likely contribute more to vaccine-induced protection against variant SARS-CoV-2 viruses.

DISCUSSION
In this study, we evaluated the efficacy of two vaccines, mRNA-1273 and Ad26.COV2.S, against the B.1.621 variant of SARS-CoV-2 in three pre-clinical models: 129S2-immuno-competed mice, K18-hACE2 transgenic mice, and Syrian hamsters. The mRNA-1273 vaccine induced high levels of neutralizing antibodies against WA1/2020 N501Y/D614G and B.1.621 viruses, and this response was associated with robust protection from an intranasal challenge. Immunization of 129S2  Our studies provide a comparison of the immunogenicity and efficacy of the Ad26. COV2.S and mRNA-1273 vaccines in 129S2 mice. The serum neutralizing antibody titer was similar after one dose of Ad26.COV2.S or two doses of mRNA-1273 in129S2 mice, and this was true for both the high-dose (5 mg versus 10 10 virus particles) and the low-dose (0.25 mg versus 10 8 virus particles) vaccine regimens. Despite the similarity in neutralization titer, the mRNA-1273 vaccine more effectively reduced viral load in the lungs of WA1/2020 N501Y/D614G-or B.1.621-challenged animals than the Ad26.COV2.S vaccine. In mice immunized with a low dose of the mRNA-1273 vaccine, approximately 20% of the animals showed evidence of breakthrough infections after challenge ( Figure 3G). In contrast, 75% of the animals that received a low dose of Ad26.COV2.S and were challenged with WA1/2020 N501Y/D614G or B.1.621 showed virus breakthrough despite relatively equivalent levels of serum neutralizing antibody titers at the time of challenge ( Figure 1F). One explanation for this difference could that the mRNA vaccine requires two doses, while the Ad26.COV2.S vaccine was only given once. Another possibility is the time between the last immunization and the virus challenge, which is 41 and 115 days for the mRNA-1273-and Ad26.COV2.S-immunized animals, respectively. Differences in the glycosylation pattern or the immunoglobulin G (IgG) subclass of antibodies between anti-S antibodies induced by mRNA-1273 and Ad26.COV2.S induced could contribute to differences in protection, as seen in non-human primates and humans. 25,26 It is also possible that mRNA-1273 vaccine induced a better anamnestic B or T cell response in 129S2 mice compared with the Ad26.COV2.S vaccine at the time of virus challenge.
The antibody response after vaccination varied between the mouse and hamster models. Immunization of mice with 5 mg of mRNA-1273 or 10 10 Ad26.COV2.S induced serum neutralizing antibody responses with a GMT of >10,000 in both 129S2 and K18-hACE2 mice. In contrast, in Syrian hamsters, the GMT against WA1/2020 was $10-fold lower than in mice yet still severalfold higher than that observed in humans vaccinated with one dose of Ad26.COV2.S. [27][28][29][30] The reason for this difference in vaccine response between mice and hamsters remains unknown. It is possible that the hamster immune response targets different epitopes on the S protein. Alternatively, the S protein contains fewer T cell epitopes for hamsters compared with mouse, although that seems unlikely given the size of the antigen. The Syrian hamster also may be more tolerogenic, perhaps due to its complex microbiome as opposed to the SPF microbiome of 129S2 and K18 TG mice. Mice that received the microbiome from pet stores or from field mice had blunted vaccine responses compared with laboratory-housed mice in pathogen-free facilities. 31 It is also possible that the differences in antibody response after immunization are impacted by differential immunological maturity between 7-and 8-week-old mice and 5-to 6-week-old hamsters. Additional studies are required to elucidate the causes for this difference, but in general, the magnitude of the antibody response in Syrian hamsters is more similar to the human antibody response after vaccination.
The B.1.621 (Mu) variant of SARS-CoV-2 has R346K and E484K mutations in the receptor-binding domain of the S protein and is believed to be more resistant to virus neutralization by serum antibodies compared with the historical SARS-CoV-2 virus. In sera from vaccinated or infected individuals, the fold difference in neutralization between the D614G (B.1) variant of SARS-CoV-2 and B.1.621 was between 2-and 12-fold. [12][13][14][15][16]32 In Syrian hamsters, we observed a $3-fold decrease (p < 0.001) in serum neutralization titers between the WA1/2020 and B.1.621 virus. In K18-hACE2 mice immunized with mRNA-1273, the difference in serum neutralization titers between WA1/2020 N501Y/D614G and B.1.621 was 1.3-to 2.3-fold. In contrast, no difference in neutralization titers between WA1/2020 N501Y/ D614G and B.1.621 was observed in 129S2 mice immunized with mRNA-1273 or Ad26.COV2.S, suggesting mouse-strain and species-specific differences in the antibody response. Another reason may be the N501Y mutation in WA1/2020 N501Y/D614G, which has been shown to reduce the neutralization titer of serum and certain mAbs compared with WA1/2020 and WA1/2020 D614G, respectively. 6,33 We used the WA1/2020 N501Y/D614G virus because the N501Y mutation was required for virus infection in immunocompetent 129S2 mice lacking hACE2 expression. Finally, it is possible that the insertion of a threonine at positions 144-145 in our particular B.1.621 isolate reduced the resistance to serum neutralizing antibodies.
We observed no difference in efficacy of the mRNA and adenoviral-vectored vaccine to protect against B.1.621 in three different animal models. Near full protection, defined by undetectable levels of SARS-CoV-2 viral RNA and infectious virus, plus the absence of immunopathology in the vaccinated animals, was observed against both B.1.621 and the control WA1/2020 virus in mice immunized with mRNA-1273. While immunization with lower doses of Ad26.COV2.S offered only partial protection, we did not observe a difference in virus titer or frequency of breakthrough infection between B.1.621 and control virus. This suggests that a 2-to 3-fold reduction in serum neutralization titers has limited impact on mRNA and adenoviral vaccine protection against the variant B.1.621 virus.

Limitations of the study
We note several limitations of our study. (1) We did not evaluate the effects of the vaccine on the transmission of SARS-CoV-2 in Syrian hamsters, which may be an important measure of vaccine protection. (2) We used lower doses of vaccine to mimic suboptimal and possibly waning immunity. Studies that directly compare the quality of a waning immune response with that of a low-dose-vaccine-induced immune response are needed. (3) The challenge dose of SARS-CoV-2 used in our hamster model (10 3 PFU) is several orders of magnitude higher that the minimal infectious dose (5 PFU). 34 While this creates a very robust virus challenge model, it could underestimate the protective effects of vaccines. (4) We did not establish correlates of immune protection. We noted that lower serum antibody neutralization titers were associated with high viral loads and infectious virus titers in the Ad26.COV2.S-immunized animals, as well as breakthrough infections in some of the mRNA-1273 vaccinated mice, albeit this did not explain all breakthrough infections. A more detailed analysis of T cell and non-neutralizing antibody responses coupled with even lower vaccine doses may be needed to fully establish a correlate of protection against breakthrough infection.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Materials availability
All requests for resources and reagents should be directed to the lead contact. This includes viruses, vaccines, and primer-probe sets. All reagents will be made available on request after completion of a Materials Transfer Agreement.
Data and code availability All data supporting the findings of this study are available within the paper and are available from the corresponding author upon request. This paper does not include original code. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS Cells and viruses
Vero cells expressing human ACE2 and TMPRSS2 (Vero-hACE2-hTMPRSS2, 6 has been published previously 6 and was propagated on Vero-hTMPRSS2 cells.
The virus stocks were subjected to next-generation sequencing, and the S protein sequences were identical to the original isolates. The infectious virus titer was determined by plaque or focus-forming assay on Vero-hACE2-hTMPRSS2 or Vero-hTMPRSS2 cells.
Pre-clinical vaccine mRNA and lipid nanoparticle production process A sequence-optimized mRNA encoding prefusion-stabilized Wuhan-Hu-1 (mRNA-1273) SARS-CoV-2 S-2P protein was synthesized in vitro using an optimized T7 RNA polymerase-mediated transcription reaction with complete replacement of uridine by N1m-pseudouridine. 37 The reaction included a DNA template containing the immunogen open-reading frame flanked by 5 0 untranslated region (UTR) and 3 0 UTR sequences and was terminated by an encoded polyA tail. After transcription, the cap-1 structure was added to the 5 0 end using the vaccinia virus capping enzyme (New England Biolabs) and vaccinia virus 2 0 -O-methyltransferase (New England Biolabs). The mRNA was purified by oligo-dT affinity purification, buffer exchanged by tangential flow filtration into sodium acetate, pH 5.0, sterile filtered, and kept frozen at À20 C until further use.
The mRNA was encapsulated in a lipid nanoparticle through a modified ethanoldrop nanoprecipitation process described previously. 38 Ionizable, structural, helper, and polyethylene glycol lipids were briefly mixed with mRNA in an acetate buffer, pH 5.0, at a ratio of 2.5:1 (lipid:mRNA). The mixture was neutralized with Tris-HCl, pH 7.5, sucrose was added as a cryoprotectant, and the final solution was sterilefiltered. Vials were filled with formulated lipid nanonparticle and stored frozen at À20 C until further use. The vaccine product underwent analytical characterization, ll e2 Med 3, 309-324.e1-e6, May 13, 2022 which included the determination of particle size and polydispersity, encapsulation, mRNA purity, double-stranded RNA content, osmolality, pH, endotoxin, and bioburden, and the material was deemed acceptable for in vivo study.

Mouse experiments
Animal studies were carried out in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocols were approved by the Institutional Animal Care and Use Committee at the Washington University School of Medicine (assurance number A3381-01). Seven-to-nine week old male 129S2 (strain: 129S2/SvPasCrl, Cat # 287) or female K18-hACE2 transgenic mice (strain: 2B6.Cg-Tg(K18-ACE2)2Prlmn/ J, Cat # 34860) were obtained from Charles River Laboratories and Jackson Laboratories, respectively and housed at Washington University. Animals were housed in groups and fed standard chow diet.
Some of the animals were vaccinated via intramuscular (IM) route with 10 8 , 10 9 , or 10 10 viral particles of fresh or freeze-thawed Ad26.COV2.S in 100 mL of phosphate buffered saline (PBS). The freeze-thawed vaccine was stored at -80 C prior to thawing at room temperature. No difference in serum antibody responses were detected between the fresh and freeze-thawed Ad26.COV2.S vaccine (Figures S1A and S1B). Control animals for the adenoviral vaccine received PBS alone. Twenty-one days and 115 days later, serum samples were obtained, and used for ELISA and virus neutralization assays. This later time point was selected to mimic potential waning immunity. Separately, 129S2 mice and K18-hACE2 mice were immunized and boosted with 0.25 or 5 mg of mRNA-1273 or a control mRNA (mRNA-control) vaccine at three week intervals. Twenty-one days after the second immunization, serum was obtained and used for virus neutralization assays. These two mouse models were selected because of their different characteristics: (a) K18-hACE2 mice are among the most susceptible model for SARS-CoV-2 infection and accounts for the SARS-CoV-2 spike protein-hACE2 entry pathways and the impact of vaccine induced antibodies on this interaction; and (b) 129 mice were used because in contrast to BALB/c and C57BL/6 mice, in 6 to 8 week-old animals, they lose weight upon challenge with 10 5 PFU of SARS-CoV-2.
Following transfer to the enhanced Biosafety level 3 laboratory, the animals were challenged via intranasal route with 10 3 or 10 5 PFU of the SARS-CoV-2 N501Y/ D614G or B.1.621 variant. This differences in challenge dose reflects the differences in susceptibility to SARS-CoV-2 between the 129S2 (10 5 PFU) and K18-hACE2 ll transgenic mice (10 3 PFU). Animal weights were measured daily for the duration of the experiment. At different time points after challenge, the animals were necropsied, and their lungs were collected for virological and histological analysis. These time points for necropsies and tissue collections were based on optimized ones that were established in the Boon and Diamond laboratories. 19,21,22,24,41,42 These time points coincided with maximal infection, pathological findings, and weight loss. The left lobe was homogenized in 1.0 mL of Dulbecco's Modified Eagle Medium (DMEM), clarified by centrifugation (1,000 x g for 5 min) and used for viral titer analysis by quantitative RT-PCR (RT-qPCR) using primers and probes targeting the N gene, and by plaque assay. A nasal wash also was collected, by inoculating 1.0 mL of PBS with 0.1% bovine serum albumin into one nostril and collecting the wash from the other nostril (Figures 1 and 4). Alternatively, 0.5 mL of PBS with 0.1% bovine serum albumin was flushed through the nasal cavity after dissecting off the lower jaw (Figures 2 and 3). The nasal wash was clarified by centrifugation (2,000 x g for 10 min) and used for viral titer analysis by RT-qPCR using primers and probes targeting the N gene, and by plaque assay.

Hamster experiments
Animal studies were carried out in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocols were approved by the Institutional Animal Care and Use Committee at the Washington University School of Medicine (assurance number A3381-01). Five-week old male hamsters were obtained from Charles River Laboratories and housed at Washington University. Five days after arrival, the animals were immunized via intramuscular injection with 10 8 of 10 10 viral particles of freeze-thawed Ad26.COV2.S in 100 mL of PBS. Control animals received PBS alone. Serum samples were obtained 21 days later, and six weeks after that the animals were transferred to the enhanced Biosafety level 3 laboratory. This time point was selected to assess responses in the setting of possible waning immunity. One day later, the animals were challenged via intranasal route with 10 3 PFU of WA1/2020 or B.1.621 variant. Syrian hamsters are highly susceptible to SARS-CoV-2 with an estimated infectious dose of 5 virus particles. 34 Therefore, we chose a challenge dose that is approximately 100-fold higher to evaluate vaccine efficacy in this model. Animal weights were measured daily for the duration of the experiment. Four days after challenge, the animals were necropsied, and their lungs were collected for virological and histological analysis. The left lobe was homogenized in 1.0 mL of DMEM, clarified by centrifugation (1,000 x g for 5 min) and used for viral titer analysis by quantitative RT-PCR using primers and probes targeting the N gene, and by plaque assay. A nasal wash was also collected, by inoculating 1.0 mL of PBS with 0.1% bovine serum albumin into one nostril and collecting the wash from the other nostril. The nasal wash was clarified by centrifugation (2,000 x g for 10 min) and used for viral titer analysis by quantitative RT-PCR using primers and probes targeting the N gene, and by plaque assay. SARS2-57, and SARS2-71 43 anti-S protein antibodies and HRP-conjugated goat anti-mouse IgG (Sigma Cat # A8924) in PBS supplemented with 0.1% saponin and 0.1% bovine serum albumin. SARS-CoV-2-infected cell foci were visualized using TrueBlue peroxidase substrate (KPL) and quantitated on an ImmunoSpot microanalyzer (Cellular Technologies).

Virus titration assays
Plaque assays were performed on Vero-hACE2-hTRMPSS2 cells in 24-well plates. Lung tissue homogenates or nasal washes were diluted serially by 10-fold, starting at 1:10, in cell infection medium (DMEM + 2% FBS + 100 U/mL of penicillin-streptomycin). Two hundred and fifty microliters of the diluted virus were added to a single well per dilution per sample. After 1 h at 37 C, the inoculum was aspirated, the cells were washed with PBS, and a 1% methylcellulose overlay in MEM supplemented with 2% FBS was added. Seventy-two hours after virus inoculation, the cells were fixed with 4% formalin, and the monolayer was stained with crystal violet (0.5% w/v in 25% methanol in water) for 1 h at 20 C. The number of plaques were counted and used to calculate the plaque forming units/mL (PFU/mL).
To quantify viral load in lung tissue homogenates and nasal washes, RNA was extracted from 100 mL samples using E.Z.N.A.â Total RNA Kit I (Omega) and eluted with 50 mL of water. Four microliters RNA was used for real-time RT-qPCR to detect and quantify N gene of SARS-CoV-2 using TaqManä RNA-to-CT 1-Step Kit (Thermo Fisher Scientific) as described 44 using the following primers and probes: Forward: GACCCCAAAATCAGCGAAAT; Reverse: TCTGGTTACTGCCAGTTGAAT CTG; Probe: ACCCCGCATTACGTTTGGTGGACC; 5 0 Dye/3 0 Quencher: 6-FAM/ ZEN/IBFQ. Viral RNA was expressed as N gene copy numbers per mg for lung tissue homogenates or mL for nasal swabs and nasal washes, based on a standard included in the assay, which was created via in vitro transcription of a synthetic DNA molecule containing the target region of the N gene.

Histology
The lungs from SARS-CoV-2 infected and control mice were inflated and fixed in 10% formalin for seven days. Lungs were embedded in paraffin and sectioned before hematoxylin and eosin staining. Lung slides were scanned using the Hamamatsu NanoZoomer slide scanning system and head sections were imaged using the Zeiss AxioImager Z2 system.

ELISA
Ninety-six-well microtiter plates (Nunc MaxiSorp; ThermoFisher Scientific) were coated with 100 mL of recombinant SARS-CoV-2 S protein (Wuhan strain) at a concentration of 1 mg/mL in PBS (Gibco) at 4 C overnight; negative control wells were coated with 1 mg/mL of BSA (Sigma). Plates were blocked for 1.5 h at room temperature with 280 mL of blocking solution (PBS supplemented with 0.05% Tween-20 (Sigma) and 10% FBS (Corning)). Serum from mice and hamsters were diluted serially in blocking solution, starting at 1:100 dilution and incubated for 1.5 h at room temperature. The plates were washed three times with T-PBS (1X PBS supplemented with 0.05% Tween-20), and 100 mL of goat anti-mouse IgG (Southern Biotech Cat #1030-05) diluted 1:2,000 in blocking solution or 100 mL of HRP-conjugated antihamster IgG(H+L) antibody (Southern Biotech Cat. #6061-05) diluted 1:500 in blocking solution, was added to all wells and incubated for 1 h at room temperature. Plates were washed 3 times with T-PBS and 3 times with 1X PBS, and 100 mL of 1-step Ultra TMB-ELISA substrate solution (Thermo Fisher Scientific) was added to ll all wells. The reaction was stopped after 5 min using 100 mL of 1M HCl, and the plates were analyzed at a wavelength of 490 nm using a microtiter plate reader (BioTek).

QUANTIFICATION AND STATISTICAL ANALYSIS
Statistical significance was assigned when P values were < 0.05 using GraphPad Prism version 9.3. Tests, number of animals, median and geometric mean values, and statistical comparison groups are indicated in the Figure legends. Analysis of weight change was determined by two-way ANOVA. Changes in infectious virus titer, viral RNA levels, or serum antibody responses were compared to unvaccinated or mRNA-control immunized animals and were analyzed by one-way ANOVA with a multiple comparisons correction, unpaired t-test, or Mann-Whitney test, dependent on the number of comparison and the distribution of the data.


A New Strategy for Multitarget Drug Discovery/Repositioning Through the Identification of Similar 3D Amino Acid Patterns Among Proteins Structures: The Case of Tafluprost and its Effects on Cardiac Ion Channels

The identification of similar three-dimensional (3D) amino acid patterns among different proteins might be helpful to explain the polypharmacological profile of many currently used drugs. Also, it would be a reasonable first step for the design of novel multitarget compounds. Most of the current computational tools employed for this aim are limited to the comparisons among known binding sites, and do not consider several additional important 3D patterns such as allosteric sites or other conserved motifs. In the present work, we introduce Geomfinder2.0, which is a new and improved version of our previously described algorithm for the deep exploration and discovery of similar and druggable 3D patterns. As compared with the original version, substantial improvements that have been incorporated to our software allow: (i) to compare quaternary structures, (ii) to deal with a list of pairs of structures, (iii) to know how druggable is the zone where similar 3D patterns are detected and (iv) to significantly reduce the execution time. Thus, the new algorithm achieves up to 353x speedup as compared to the previous sequential version, allowing the exploration of a significant number of quaternary structures in a reasonable time. In order to illustrate the potential of the updated Geomfinder version, we show a case of use in which similar 3D patterns were detected in the cardiac ions channels NaV1.5 and TASK-1. These channels are quite different in terms of structure, sequence and function and both have been regarded as important targets for drugs aimed at treating atrial fibrillation. Finally, we describe the in vitro effects of tafluprost (a drug currently used to treat glaucoma, which was identified as a novel putative ligand of NaV1.5 and TASK-1) upon both ion channels’ activity and discuss its possible repositioning as a novel antiarrhythmic drug.

INTRODUCTION
Although most novel drugs are still developed using the "magic bullet" paradigm, which involves highly selective profiles, chemical compounds are naturally promiscuous in practice. Indeed, most therapeutically beneficial agents interact with more than one molecular target (Feldmann and Bajorath, 2020). Interestingly, this promiscuity is now, in some cases, considered an advantageous feature and is proactively pursued. Thus, many drug development initiatives are focused on the design of multitarget compounds with a polypharmacological profile that may improve drug-based treatments' efficacy and/or safety (Ramsay et al., 2018;Proschak et al., 2019). Unfortunately, the design of drugs with multiple activities on a selected handful of different protein structures remains a significant experimental and computational challenge (Konc, 2019). Recent reports have proposed several strategies to identify multitarget drugs, such as clinical observations and target combinations based on phenotypic screening. Bioinformatics is also a useful tool to address these challenges through molecular modelling techniques for detecting similar targets, machine learning methods to find disease-related targets, target-fishing using molecular docking, ligand-based pharmacophore searching, virtual screening simulations, and the search of binding sites similarities (Ma et al., 2010;Ren et al., 2021;Stępnicki et al., 2021). Finding compounds with multitarget action on related proteins which share a similar function, folding or binding sites, is currently an accessible task. Unfortunately, complex diseases often comprise a wide range of evolutionary distant and structurally different proteins where current methods are not entirely precise. For example, neuropsychiatric, cardiac or autoimmune disorders (among others), including atrial fibrillation or major depression, are complex diseases that often encompass dysfunctions in a wide range of types of proteins such as ion channels, enzymes, transporters, globular proteins, etc. (Bolognesi, 2019;Konc, 2019). Thus, the identification of conserved/similar sites (or more broadly three-dimensional (3D) patterns, defined as a local structural arrangement of amino acids) among a set of proteins (related or not between them), can be useful for the rational design/ repositioning of polypharmacological drugs (Konc, 2019;Adasme et al., 2020;Li et al., 2021). In this context, tools such as G-LoSA (Lee and Im, 2016), Geomfinder (Núñez-Vivanco et al., 2016), 3D-PP (Valdés-Jiménez et al., 2019) and others (Ehrt et al., 2016;Ehrt et al., 2018), which work regardless of information about known ligands, binding sites, sequence similarity, or structural folding of the proteins, improve the chances of finding similar 3D patterns among very different or unrelated targets. Although some of these similarities may appear by chance, others might represent distant evolutionary relationships and correspond, for instance, to secondary binding sites. These sites have recently gained attention for the rational design of polypharmacological allosteric modulators (Abdel-Magid, 2015;Meysman et al., 2015;Reyes-Parada and Iturriaga-Vasquez, 2016;Wakefield at al., 2019). Indeed, by using our in-house algorithm Geomfinder (Núñez-Vivanco et al., 2016), we have reported the finding of some similar 3D patterns among very different protein structures, which cannot be observed through other structural tools or with sequencebased methods. Noteworthy, in the original version of Geomfinder, the residues of each detected 3D pattern could only be part of one specific protein chain. However, it is well known that several important 3D patterns (e.g. binding sites, catalytic sites) are located either at the interface between different subunits of a single target (Lee et al., 2017) or at the oligomer interface in multimeric proteins (Baek et al., 2017).
In the present work, we introduce Geomfinder2.0, which is a new and improved version of our previous tool for the deep exploration and discovery of similar and druggable 3D patterns. Thus, the possibility of exploring 3D patterns formed by different chains of a quaternary protein structure is one of the major novel features of this version. Also, the updated algorithm includes a function that predicts the zone where similar 3D patterns may be druggable (Le Guilloux et al., 2009). Thus, a 3D pattern with a high level of druggability found in several protein structures might be used as the input for the design of multitarget compounds using approaches such as Pocket-Based Drug Design (Zheng et al., 2013). From a computational perspective, the currently available version of Geomfinder has been fully migrated from Python 2.7 to C++ language and parallelized using the shared memory programming model OpenMP (Dagum and Menon, 1998). These improvements allowed a speedup of up to 353x. Furthermore, as a functional enhancement, Geomfinder can now compute several pairs of comparisons simultaneously.
In addition, in order to illustrate the potential of the updated Geomfinder version, we show a case of use in which similar 3D patterns were detected in two cardiac ion channels, specifically NaV1.5 and TASK-1, which are selective for sodium or potassium, respectively. These channels have been regarded as important targets for drugs aimed at treating atrial fibrillation (Sossalla et al., 2010;Wiedmann et al., 2021). Finally, we describe the in vitro effects of tafluprost (a putative ligand identified after a receptor-based drug search approach) upon both ion channels' activity and discuss its possible repositioning as a novel antiarrhythmic drug.

Software Improvements
The new version of Geomfinder includes all geometrical characteristics and processes described previously (Núñez-Vivanco et al., 2016) and incorporates new features and substantial improvements. In the Supplementary Figure S1, the implemented architecture and essential components and services of Geomfinder are shown. Remarkably, the input can be a list of pairs of structures (always as PDB files) in this new version. This new feature makes it possible to find similar 3D patterns simultaneously in several pairs of protein structures submitted by the user (Figure 1). Another essential feature of the new version is the possibility of searching for 3D patterns in the interface between two or more protein chains/subunits. To   March 2022 | Volume 13 | Article 855792 3 this end, in the submission process, the user can select one particular chain of each or concatenate either all or some of the structure's chains ( Figure 2). Also, it is now possible to search 3D patterns only in those zones of the proteins where the cavities detected achieve the user-defined druggability score threshold, which is interpreted as the probability of the cavity to bind a drug and alter its normal activity (Pérot et al., 2010;Schmidtke and Barril, 2010). This preprocessing step is an optional parameter ( Figure 3) calculated by the Fpocket algorithm (Le Guilloux et al., 2009), considering features such as the size, the hydrophobicity, and the normalized polarity of the residues lining the cavity. Technically, several sections of the algorithm were also optimized. The original version of Geomfinder (PythonThreading) was developed in Python 2.7 using threads as a parallel strategy. Searching for a better performance, three new versions were built: PythonMultiprocessing using Python 2.7 with multiprocessing, C++/Pthreads using C++ and POSIX threads, and C++/OpenMP using C++ with OpenMP annotations. These versions were analyzed (using a 32-CPU machine with hyperthreading activated) against a sequential version of Geomfinder implemented with benchmarking purposes.

Detection of Similar 3D Amino Acid Patterns Using Geomfinder
Searching for similar 3D amino acid patterns begins with the creation of a virtual grid of coordinates on each protein structure. Then, using all geometrical centers of all side chains of the residues, Geomfinder performs a residues grouping step depending on the distance between each virtual coordinate. With this, hundreds of 3D amino acid patterns are defined in each protein structure. After that, for each 3D pattern, four descriptors are measured: a) The list of distances between the geometric centers of the side chains of all the residues forming the 3D pattern; b) the sum of the short and medium-range of non-bonded energy of each residue forming the 3D pattern; c) the list of the residues forming a 3D pattern; and d) the list of the distances constituting the shortest pathway necessary to go over all the residues lining the 3D pattern. Finally, all pairs of 3D patterns identified in the two tested proteins are compared using an all-versus-all approach. Thus, at the end of the analysis, each pair of the 3D pattern has a final similarity score named GScore. The GScore is defined as a combination of the similarities (S) of the four descriptors: which is calculated as the relative changes on each pair of the 3D pattern as follows: SDist, SNbE, STsp, and SSc, are the partial scores of similarity of the distances, the non-bonded energies, the perimeter, and the sequence components, between any two 3D patterns.

Receptophore Determination
As we have previously proposed (Núñez-Vivanco et al., 2018), a "receptophore" can be defined -by analogy with the pharmacophore concept-as a 3D ensemble (present in two or more receptors), of molecular, steric, and electronic features that ensure the optimal molecular interactions with a common promiscuous ligand. Therefore, here we describe how we determine it from the local similarities identified with Geomfinder. The method consists basically of the structural alignment of the similar 3D patterns identified between the protein structures. This process was performed using the external computational methods PocketAlign and MultiBind (Shulman-Peleg et al., 2008;Yeturu and Chandra, 2011). This approach finds the best match of physicochemical properties among the residues forming each site. PocketAlign carries out multiple alignments between the similar 3D patterns detected to recognize similar amino acids matches. Then, numerous structural rearrangements of superimposed binding sites are applied to find the best structural fit. Briefly, this method consists of two main processes: a) the preprocessing of the features of each 3D pattern and hashing them into a table; and b) the recognition of the similar features in the objects of the hash table. In the preprocessing, each amino acid is denoted by pseudocenters (X, Y, and Z coordinates), which provides a unique physicochemical property to the binding site: hydrogenbond donor, hydrogen-bond acceptor, mixed donor/acceptor, hydrophobic aliphatic or aromatic contacts. Finally, MultiBind performs a combination of multiple superimposed binding site conformations to find common patterns. Because 3D patterns that originate the receptophore are not necessarily identical, different residues might be partially aligned if the overall structural alignment score turns out to be maximized with that arrangement. After that, to generate the final receptophore, all equivalent amino acids between both 3D patterns (same physicochemical group: polar, non-polar, positively or negatively charged) appearing aligned or superimposed, are manually merged in the resultant PDB file. In contrast, all non-equivalent amino acids are preserved in the final receptophore.

Receptor/Pocket Based Virtual Screening
The previously defined receptophore was used to perform a virtual screening analysis through the computation tool e-LEA3D (Douguet, 2010). This program can create new molecules with a fragment-based approach or evaluate a user-defined data set of compounds. In our case, we used the "FDA-approved" data set to determine if some of these compounds exhibit affinity for the receptophore. Thus, each molecule's fitness was evaluated via a function that inputs the molecular structure and returns a numeric score. The evaluation can integrate a selected number of molecular properties and/or a protein-ligand docking score calculated by the program PLANTS (Korb et al., 2009). After that, a list of FDA-approved drugs, ranked by their affinity for the receptophore, was obtained and analyzed for further experiments.

Protein Structures
For the case of use, we employed the crystallographic protein structures (obtained from the Protein Data Bank) of the sodium channel NaV1.5 (PDBid: 6UZ3, resolution 3.50 Å) and the twopore domain potassium ion channel TASK-1 (PDBid: 6RV3, resolution 2.90 Å). Only chain A was selected in the case of NaV1.5, whereas for TASK-1, chains A and B were concatenated. All details of input parameters and results can be found at https:// geomfinder2.appsbio.utalca.cl/result/11598984329232/

Drugs and IC50 Values
Tafluprost was purchased from Merck. All other reagents used were of analytical grade. Tafluprost was dissolved in DMSO and added to the external solution just before the recordings. The IC50 (or EC50) were determined from Hill plots using up to five concentrations for each construct and are expressed as mean ± SD coming from the different replicates measurements (n = 3-9 replicates).

Oocyte Preparation and cRNA Injection
All procedures performed in this study involving animals were carried out in accordance with the EU Directive 2010/63/EU for animal experiments. The work with Xenopus laevis at the University of Marburg with all experimental protocols were approved by the Regierungspräsidium Gießen, Germany (V54-19c20 15 h 02 MR 20/28 Nr.A 4/2013).
Oocytes were obtained from anesthetized Xenopus laevis frogs and incubated in OR2 solution containing in mM: 82.5 NaCl, 2 KCl, 1 MgCl2, 5 HEPES (pH 7.5) substituted with 2 mg/ml collagenase II (Sigma) to remove residual connective tissue. Then the oocytes were stored at 18°C in ND96 supplemented with 50 mg/L gentamycin, 274 mg/L sodium pyruvate, and 88 mg/L theophylline.
For both, TASK-1 or NaV1.5 channel measurements, recording solution ND96 contained in mM: 96 NaCl, 2 KCl, 1.8CaCl2, 1 MgCl2, 5 HEPES (pH 7.5). In the case of TASK-1 channel, block was analyzed with voltage steps from a holding potential of −80 mV. A first test pulse to 0 mV of 1 s duration was followed by a repolarizing step to −80 mV for 1 s directly followed by another 1 s test pulse to +40 mV. The sweep time interval was 10 s. For NaV1.5 channel, block was analyzed with voltage steps from a holding potential of −120 mV. A first depolarizing pulse to −10 mV of 20 ms duration, then, after holding at −40 mV for 4s, a 20 ms step at −120 mV was carried out. The sweep time interval was 10 s.
Tafluprost at different concentrations was evaluated on channel currents. Stability in recordings was monitored prior to the addition of compounds, which were removed from the bath to show recovery.

Software Optimization
All the new versions of Geomfinder, implemented with benchmarking purposes, incorporated several optimizations such as more efficient data structures and compilation flags for the machines utilized on the webserver. As denoted by the green line in Figure 4, the version C++/OpenMP always showed a better performance than the other implementations, achieving a maximal 353x speedup compared to the original Sequential version (Supplementary Figure S6). On this basis, the C++/ OpenMP version was selected for the implementation of the new server of Geomfinder.

Case of Use
It has been shown that some local anesthetics are multichannel blocking drugs, which interact with cardiac ion channels such as NaV1.5 and TASK-1 (Tikhonov and Zhorov, 2017;Rinné et al., 2019). Considering that this polypharmacological profile likely underlies their antiarrhythmic effect, we searched similar and druggable 3D patterns between these channels. Although NaV1.5 and TASK-1 have different sequences, structures, and topologies ( Figure 5), Geomfinder required a few seconds to detect several similar 3D patterns (Supplementary Figure S2).
Interestingly, one of the patterns found contains two key residues of the local anesthetics binding site in NaV1.5 (Phe1762, Ile1468; (Nguyen et al., 2019)) and one key residue in TASK-1 (Phe238; (Rinné et al., 2019)) located at the fenestrations. Considering that these 3D patterns are similar (GScore = 66.5%) and are located in zones with high druggability values (Phe1762 and Ile1468 are located in pockets with 99.9% druggability in NaV1.5; and Phe238 in a pocket with 100% druggability in TASK-1; result #25 in Supplementary Figure S3), it is enticing to state that the use of these 3D patterns is a promising starting point to understand the polypharmacological profile of local anesthetics and to design/search new multitarget compounds aimed to these cardiac ion channels. The fact that the residues Phe1762 in NaV1.5 and Phe238 in TASK-1 were found by Geomfinder in similar tridimensional orientations into each 3D pattern (Figure 6), supports the idea that the aromatic ring of local anesthetics interacts with the ion channels establishing a π−π interaction (González et al., 2001). FIGURE 4 | Performance of new Geomfinder implementations. The original (PythonThreading) and three new versions, including parallel programming paradigm and code optimization, were compared against a sequential Python version (Sequential) developed with benchmarking purposes. These versions were built using Python 2.7 with multiprocessing (PythonMultiprocessing), C++ with POSIX threads (C++/Pthreads), and C++ with OpenMP annotations (C++/OpenMP). As is denoted by the green line in the graph, the version C++/OpenMP always showed better performance than the other implementations, reaching their best acceleration using 54 threads: 353x speedup compared to Sequential version. Thus, the version C++/OpenMP was selected for the implementation of the new server of Geomfinder. The experimental setup used in the comparisons is a 32 Intel Xeon CPU E5-2683 (2.10 GHz) SMP system with hyperthreading enabled (64 virtual cores/threads), 252 GB RAM, 40 MB Intel Smart Cache.
FIGURE 5 | General overview of structure (A), topology (B) and sequence (inset) similarity of NaV1.5 and TASK-1. NaV1.5 corresponds to the PDBid:6UZ3 and TASK-1 to the PDBid:6RV3. NaV1.5 is a monomer with four domains, each containing a pore sequence and TASK-1 is a dimer with four transmembrane segments and two pore sequences each monomer. In order to show how these results could be useful for the design/search of new multitarget compounds, we initially performed docking molecular simulations of bupivacaine (a common local anesthetic) on NaV1.5 and TASK-1 using the coordinates of the similar pockets detected with Geomfinder as grid centers. As shown in Supplementary Figures S4, S5, the aromatic ring of the bupivacaine seems to be establishing π-π interactions with Phe238 in TASK-1 and with Phe1762 in NaV1.5. After selecting the best conformers of bupivacaine (those with the lowest estimation of free energy of binding) in both proteins, we extracted the residues located at 4 Å of the ligands, and constructed a common binding site for TASK-1 and NaV1.5, as we have previously described (Möller-Acuña et al., 2015;Núñez-Vivanco et al., 2018). Briefly, this common binding site was constructed using external computational tools such as PocketAlign (Yeturu andChandra, 2011) andMultiBind (Shulman-Peleg et al., 2008), which perform several structural alignments and tridimensional rearrangements, looking for the best match of the physicochemical properties of the selected binding sites. As shown in Figure 7, it was possible to find a tridimensional fit for adjusting the physicochemical properties and the structural alignment of the residues of both binding sites, FIGURE 6 | Similar 3D pattern found in NaV1.5 and TASK-1. In both sites, the residues PHE1762 and PHE238 show similar orientation. Frontiers in Pharmacology | www.frontiersin.org March 2022 | Volume 13 | Article 855792 7 which we call the common "receptophore" (Núñez-Vivanco et al., 2018). After that, we used this "receptophore" to perform a structure-based virtual screening with the e-LEA3D tool (Douguet, 2010). This software requires as input a PDB file of one pocket structure, which can then be used to make either a de novo drug design or a virtual screening into a collection of known compounds, which in both cases should lead to find molecules with high potential affinity for the pocket. In our case, we merged both aligned binding sites into one unique PDB file and after run e-LEA3D in the virtual screening mode (using the "FDA approved drugs" data set available at the e-LEA3D web server), a list of molecules ranked by their theoretical affinities for the pocket submitted was obtained. The top ranked molecule was tafluprost (Drugbank_ID: DB08819), a prostaglandin F 2α (PGF 2α )-type agonist currently used as a treatment for glaucoma and ocular hypertension (Papadia et al., 2011;Klimko and Sharif, 2019), which in theory should be, a new putative multitarget ligand of NaV1.5 and TASK-1.
To test this idea, we performed new docking molecular simulations with tafluprost on the structures of NaV1.5 and TASK-1, setting the same parameters that were used for the experiments with bupivacaine. As stated in Figure 8, tafluprost showed better affinities than those obtained for bupivacaine at both protein structures. The stabilization of the protein-ligand complexes seems to be determined by several hydrophobic interactions of residues which have been previously reported as key residues for the local anesthetics (e.g. Phe238, Leu1464 and Phe1762).

Experimental Evaluation of Tafluprost
As shown in Figure 9, tafluprost potently blocked TASK-1 outward currents in a concentration-dependent manner (EC 50 = 186 ± 40 nM; Figures 9A,B). It should be noted that since the highest tafluprost concentration used in these experiments was 10 µM ( Figure 9B), here we report EC50 instead of IC50, as the hill fit determines the concentration where tafluprost reaches 50% of its maximal effect but not 50% inhibition. As we reached 82% maximal inhibition, EC50 concentration corresponds to about 40% inhibition (not 50% inhibition). On the other hand, tafluprost also blocked NaV1.5 inward currents, although with much less potency than that observed at TASK-1 (estimated IC50 of about 76 µM Figure 9B).

DISCUSSION
The search of multitarget compounds might be a difficult task, particularly when the drugs are aimed to act at receptors with highly diverse structure and function. Based on the idea that a given compound could simultaneously interact with two (or more) relevant targets if they have similar binding sites (Jalencas and Mestres, 2013;Salentin et al., 2014;Ehrt et al., 2016;Konc, 2019;Naderi et al., 2019), one reasonable approach to find promiscuous drugs under these circumstances is to look for similar binding sites at the addressed targets. In this context, Geomfinder2.0 appears as a valuable tool since it is a fast web server for the discovery of similar and druggable 3D patterns between any pair of protein structures. This new version has significantly improved its usefulness and performance as compared with the original version, with up to 353x speedup for the input data set analyzed and the available machine, and allowing to compare a list of pairs of structures. It also identifies 3D patterns formed by different protein chains and characterizes how druggable is the zone where the 3D patterns were detected. It is important to note that, beyond these functional and performance improvements, the core method for searching and comparison of 3D patterns is the same as in the original version (Núñez-Vivanco et al., 2016). The accuracy and precision of this core method has already been compared with those of computational tools such as PocketMatch and ClickTopology (Núñez-Vivanco et al., 2016). In addition, in the present work we confirmed the reliability of our algorithm with a case of use, in which the theoretical predictions were experimentally confirmed. Interestingly, when analyzing TASK-1 and NaV1.5, two cardiac channels with highly different structures and functions, Geomfinder2.0 was able to find (in a few seconds) several similar 3D patterns, a pair of which seemed remarkably attractive since they included some key residues involved in the binding of local anesthetics at both types of ion channels (Nguyen et al., 2019;Rinné et al., 2019). On this basis, we dissected these 3D patterns and constructed a common binding site which was used to search for possible novel multitarget ligands. As the development of novel polypharmacological agents can be a difficult, timeconsuming and expensive task, drug repurposing (i.e. the establishment of new indications of existing drugs), has been proposed as an efficient alternative over the de novo drug development approach. Thus, using a "FDA-approved" data set we searched for known compounds that might show an unanticipated ability to interact with TASK-1 and NaV1.5. Tafluprost, a prostaglandin analogue (PGF 2α agonist) currently used for the treatment of open-angle glaucoma (Papadia et al., 2011;Klimko and Sharif, 2019), was the most promising drug arising from this analysis. Remarkably, when experimentally tested tafluprost blocked the corresponding currents at both NaV1.5 and TASK-1, although with a higher potency for the latter. The differential activity of tafluprost upon both ion channels (as well as its higher potency as compared with the local anesthetic bupivacaine (Stoetzer et al., 2016)) roughly agrees with the theoretical binding energies predicted by the docking simulations. Noteworthy, although tafluprost was clearly less active in NaV1.5 than in TASK-1, it still shows a potency in a similar range as that shown by the well-known sodium channel blocker bupivacaine (Zhang et al., 2014). Even though this is the first time that an activity of tafluprost on cardiac channels is described, it should be noted that it had been reported that the drug is able to induce a relaxation of rabbit ciliary arteries precontracted with a high-potassium solution (Dong et al., 2008). Accordingly, this effect might be related with its potent TASK-1 blocking properties.
Both NaV1.5 and TASK-1 are attractive drug targets in particular for the development of treatments of atrial fibrillation, the most common cardiac arrhythmia (Sossalla et al., 2010;Wiedmann et al., 2021). In addition, multichannel blockers such as amiodarone or dronedarone have been used in clinical settings . In this context, the pharmacological activity predicted and demonstrated here for tafluprost, suggests that this compound could be a good candidate to evaluate its properties and repurpose it as a novel antiarrhythmic drug.
In summary, despite mounting evidence indicating that polypharmacological drugs might show better efficacy and less side effects than more selective compounds, early work in this area had already recognized that the rational search of multitarget drugs faces at least two major challenges, including a) the need to identify a combination of nodes in a biological network whose perturbation results in a desired therapeutic outcome, and b) to develop drugs whose polypharmacological profile allows those nodes to be perturbed specifically (Hopkins, 2008). Therefore, computational tools such as Geomfinder can be helpful to discover similar and druggable 3D patterns among proteins which have been tagged as targets in a multi-factorial disease, which appears as an important first step to either rationally design or -as in this case-purpose novel indications for compounds already in use.


Distinct Roles of Estrogen Receptors in the Regulation of Vitellogenin Expression in Orange-Spotted Grouper (Epinephelus coioides)

During their breeding season, estrogen induces vitellogenin (VTG) production in the liver of teleost fish through estrogen receptors (ERs) that support oocyte vitellogenesis. There are at least three ER subtypes in teleost fish, but their roles in mediating E2-induced VTG expression have yet to be ascertained. In this study, we investigated the expression of vtgs and ers in the liver of orange-spotted grouper (Epinephelus coioides). Their expression levels were significantly increased in the breeding season and were upregulated by an estradiol (E2) injection in female fish, except for the expression of erβ1. The upregulation of vtgs, erα and erβ2 by E2 was also observed in primary hepatocytes, but these stimulatory effects could be abolished by ER antagonist ICI182780 treatment. Subsequent studies showed that ERβ antagonist Cyclofenil downregulated the E2-induced expression of vtg, erα, and erβ2, while the ERβ agonist DPN simulated their expression. Knockdown of erβ2 by siRNA further confirmed that ERβ2 mediated the E2-induced expression of vtgs and erα. To reveal the mechanism of ERβ2 in the regulation of erα expression, the erα promoter was cloned, and its activity was examined in cells. E2 treatment simulated the activity of the erα promoter in the presence of ERβ2. Deletions and site-directed mutations showed that the E2 up-regulated transcriptional activity of erα occurs through a classical half-estrogen response element- (ERE) dependent pathway. This study reveals the roles of ER subtypes in VTG expression in orange-spotted grouper and provides a possible explanation for the rapid and efficient VTG production in this species during the breeding season.

Introduction
Vitellogenin (VTG), the precursor of the yolk protein in egg-bearing vertebrates, is essential for oogenesis, embryonic development, and larval survival. During the reproductive season, Vtg in the female liver is synthesized by estrogen stimulation and secreted into the bloodstream. Plasma Vtg is taken up by Vtg receptors on the surfaces of oocytes and cleaved into smaller yolk proteins by cathepsin [1]. Vitellogenesis is the primary phase of oocyte development and consumes a significant amount of nutrients and energy [2]. An insufficient vitellogenin uptake by oocytes leads to incomplete larval development and the higher mortality of eggs [3]. In female zebrafish, the knockout of vtg genes results in a dramatic drop in egg fertilization rates and offspring survival [4]. Therefore, substantial Vtg production is critical to the reproductive success of teleosts.
Other piscine ER subtypes are also expressed in the liver besides ERα [17,18]. Whether they are involved in VTG production and the regulatory relationships between ER subtypes in this process attract attention. An intraperitoneal injection of E 2 could upregulate the mRNA expression of erβ1 and erβ2 in the liver of goldfish [19]. ERβ selective agonist DPN increased VTG expression in rainbow trout hepatocytes in a dose-dependent manner [20]. A knockdown of erβ1 or erβ2 using specific siRNA in goldfish primary hepatocytes significantly reduced the E 2 -induced mRNA expression of vtg and erα [21]. Research on zebrafish embryos showed that the knockdown of erβs via specific morpholinos could suppress the estradiol induction of vtg and erα mRNA expression [22]. These studies indicate that ERβ subtypes play important roles in the vitellogenesis of fish. It is proposed that ERβ subtypes mediate estrogen signals to promote erα expression, which then sensitize the hepatocytes to further E 2 stimulation and prepare them for vitellogenesis [21]. However, a systematic study in fish revealing the interplay between ER subtypes in VTG production and the regulatory mechanism is lacking. Moreover, further studies in different fish species with different reproductive modes are required to evaluate the significance of ER subtypes in this physiological process and reveal other important aspects of ER subtypes regarding the regulation of VTG production.
Groupers are marine benthic fish, belonging to the subfamily Epinephelinae (Teleoste: Serranidae), which contains more than 160 species in 15 genera [23]. Most groupers are protogynous hermaphrodites and are of significant economic value. According to the statistics from the Food and Agriculture Organization of the United Nations (FAO), global grouper production reached about 234,828 tons in 2019. Groupers generally enter the breeding season from the end of spring to the early autumn, and this process is mainly affected by water temperature [24]. During the breeding period, females carry 70,000 to 1,000,000 ova, depending on their size [24]. Due to large-scale farming in the Asia-Pacific region, the fertilized eggs of groupers have been in short supply over the past decade. Vitellogenesis is a crucial process for oocyte development. However, there are few studies on vitellogenesis, particularly on VTG regulation in groupers. In this study, using the orange-spotted grouper (Epinephelus coioides) as an experimental model, we aimed to evaluate the functional role of ER subtypes in the regulation of VTG production and to reveal their regulatory relationships. Firstly, the expression patterns of ers and vtgs in the liver during the breeding and non-breeding seasons were examined. Secondly, the regulation of ers and vtgs expression by E 2 , ER antagonists, and ER agonists was further investigated. Finally, the regulatory relationships between ER subtypes and the molecular mechanism of ERβ2 regulation towards erα expression were studied. Our study will not only characterize the roles of ER subtypes in the regulation of VTG production, but also reveal the transcriptional regulatory mechanism of how ERβ2 regulates erα expression in fish at the molecular level.

Expression Profiles of vtgs and ers in the Liver of Orange-Spotted Grouper
To investigate the expression profiles of vtgs and ers in the liver during vitellogenesis, a histological analysis was firstly applied to examine the ovarian developmental status of each fish. The male groupers used as controls showed normal spermatogenesis in the testis ( Figure 1A). Primary oocytes were present in the ovaries of non-breeding females ( Figure 1B), while the females in the breeding season showed vitellogenesis in the ovaries with vitellogenic oocytes ( Figure 1C). The results showed that females in the breeding season had significantly higher VTG protein and mRNA levels compared to females in the non-breeding season and males ( Figure 1D-G). All three er genes were expressed in the livers of the groupers. Their mRNA levels in females were also obviously upregulated in the breeding season ( Figure 1H-J). Evidently, erα expression showed a dramatic increase in the breeding season compared to the non-breeding season. in the liver of males and females. Data are expressed as the mean ± SEM (n = 6) and analyzed by one-way ANOVA followed by Tukey's post-hoc test. Different letters above the error bars indicate statistical differences at p < 0.05.

Effects of E 2 on the Hepatic Expression of vtgs and ers in Juvenile Female Orange-Spotted Grouper
To investigate the relationship between ers and vtgs in vivo, the juvenile females received an intraperitoneal injection with a dose of 5 mg/kg E 2 . The results showed that exogenous estrogen dramatically increased the protein and mRNA levels of vtgs in the liver compared to the control group (Figure 2A-D). However, the er genes in the liver showed distinct expression patterns in response to E 2 treatment. A remarkable increase in erα mRNA levels was observed, whereas a significant decrease in erβ1 mRNA levels was found after the E 2 administration ( Figure 2E,F). The erβ2 mRNA levels were slightly but significantly elevated following the E 2 injection ( Figure 2G). Moreover, a primary hepatocyte culture system was established using the livers of the juvenile females, and the cells were incubated with different concentrations of E 2 (from 0 to 1 µM) for 24 h. The results showed that E 2 treatment significantly stimulated the protein and mRNA expression of vtgs in a dose-dependent manner ( Figure 3A-D), and evoked the dose-dependent expression of erα and erβ2 in the hepatocytes ( Figure 3E,G). However, a dose-dependent decrease in erβ1 expression was observed with increasing concentrations of E 2 treatment ( Figure 3F).

Effects of ER Antagonists or Agonist on the Expression of vtgs and ers in Primary Hepatocytes
To reveal the roles of different ERs towards the regulation of vtg expression, ER antagonists and one agonist were applied to treat the primary hepatocytes. These drugs altered the er expression and ultimately modulated the vtg expression. As shown in Figure 4, ERα and ERβ antagonist ICI182780 markedly suppressed the increasing protein and mRNA expression of vtgs induced by E 2 in the hepatocytes and abolished the stimulatory effects of E 2 on the expression of erα and erβ2 but had no effect on the reduced expression of erβ1 caused by E 2 . Cyclofenil, an Erβ-specific antagonist, significantly attenuated the stimulatory effects of E 2 on the expression of vtgs, erα, and erβ2, but could not change the erβ1 expression patterns following co-incubation with E 2 ( Figure 5).
DPN, an ERβ specific agonist, stimulated the protein and mRNA expression of vtgs effectively, and significantly upregulated the mRNA expression of three ers in hepatocytes without the presence of E 2 ( Figure 6).

Effects of erβ2 Knockdown by siRNA on the Expression of vtgs and ers in Primary Hepatocytes
To verify whether ERβ2 could regulate the expression of erα, we established a transfection system that transfers the siRNA into primary grouper hepatocytes using an electroporation method. Flow cytometry was used to evaluate the transfection efficiency in the primary hepatocytes via electroporation. The fluorescence intensity of the hepatocytes in the control group and the siRNA-NC-Cy3 (NC) group was analyzed by flow cytometry ( Figure S1A,B), and the P1 range was defined as the strong fluorescence range. As shown in Figure S1, in the P1 range, the cell count of the NC group was significantly higher than that of the control group, indicating that the siRNA was effectively transmitted into the primary hepatocytes. To test the hypothesis wherein E 2 -induced erβ2 expression would promote the erα expression and thus accelerate the vtg expression, two distinct siRNAs (Si217 and Si353) against erβ2 were separately transfected into the hepatocytes through electroporation. As shown in Figure 7, both siRNAs could significantly reduce the E 2 -induced erβ2 expression, resulting in an about 32% and 47% reduction compared to the E 2 treatment group. As expected, the expression levels of the vtgs and erα simulated by E 2 were significantly downregulated in the hepatocytes, and the erβ1 expression was not affected by erβ2 knockdown (Figure 7). The cells were transfected with si217 or si353 (5 µg /mL), and then treated with 0.01 µM E 2 for 24 h. Data are expressed as the mean ± SEM (n = 4) and analyzed by one-way ANOVA followed by Tukey's post-hoc test. Different letters above the error bars indicate statistical differences at p < 0.05.

Estrogen Stimulation on erα Promoter Activity via ERβ2
Our results indicate that ERβ2 could mediate E2 action and then regulate the erα expression, but the transcriptional regulatory mechanism is unknown. Therefore, to reveal the mechanism of E 2 in the regulation of erα expression via ERβ2, the 5 -flanking region of the erα gene was isolated. As shown in Figure S2, the cloned sequence upstream of the TSS of erα is 1434 bp in length, containing three half-estrogen response elements (ERE) and several transcription factor binding motifs, such as Activator protein 1 (Ap1), Specificity proteins 1 (Sp1), and cAMP response element binding protein (CREB). Whether the erα promoter activity was regulated by E 2 and ERβ2 was further assessed in HEK293T cells. The results showed that E 2 treatment significantly increased the erα promoter activity in the cells co-transfected with the ERβ2 ( Figure 8A). The promoter activity was significantly elevated by E 2 at 0.1 and 1 µM in the presence of ERβ2 ( Figure 8B). To identify the E 2 response region of the erα promoter, deletion analyses were further performed. The E 2 treatment significantly increased the activity of the wild-type erα promoter in the presence of ERβ2, and the deletion of the promoter to position −1057 did not eliminate the E 2 stimulatory effect ( Figure 8C). However, deleting the promoter to position −825 could abolish the E 2 -induced promoter activity ( Figure 8C), indicating that this region is responsible for mediating the E 2 action towards triggering the transcriptional activity of the erα promoter through ERβ2. As the proximal −825 bp region of the erα promoter contains a putative Sp1 binding site and two putative half-ERE sites (marked by #1, #2) ( Figure 8D), site-directed mutagenesis analyses were carried out. The mutation of the Sp1 binding site (upstream of TSS from −953 to −944 bp) or ERE#1 site (upstream of TSS from −900 to −895 bp) did not change the E 2 -induced promoter activity in the presence of ERβ2; however, the stimulatory effect was abolished when ERE#2 (upstream of TSS from −870 to −865 bp) was mutated ( Figure 8D), indicating that ERE#2 is responsible for the E 2 / ERβ2-induced expression of erα.

Discussion
The orange-spotted grouper is a seasonally spawning fish. In the breeding season, the serum estrogen levels in females are significantly elevated [25], but the pathway of estrogen in the regulation of the VTG product in the liver via ERs remains unclear in this species. In the present study, our results showed that (1) erα was abundantly expressed in the liver in the breeding season, (2) hepatic erα expression was more sensitive to E 2 treatment than erβs in vivo and in vitro, (3) and that the E 2 -induced expression of VTG and erα was almost completely inhibited after the ICI182780 treatment, suggesting that hepatic ERα is mainly responsible for mediating the action of estrogen to promote the VTG production. These results are consistent with the findings in other fishes [26].
For erβ1, its expression in females was also significantly increased in the liver in the breeding season compared to non-breeding season. However, the E 2 treatment reduced erβ1 expression in vivo and in vitro, and this inhibitory effect was not abolished by ICI182780 and Cyclofenil. These results indicate that the elevated expression of erβ1 in the breeding season may not be induced by estrogen. A study also found that hepatic erβ1 expression was markedly decreased in zebrafish after 48 h of E 2 exposure [9]. However, the hepatic erβ1 expression level varied greatly in response to E 2 in goldfish. Some studies observed an induction of hepatic erβ1 expression by E 2 in juvenile and female goldfish, while others found a downregulation of erβ1 expression in males after E 2 implantation, suggesting that Erβ1 possesses hepatic functions in teleosts, but its role in the regulation of VTG production is unclear.
Compared with other ER subtypes, erβ2 showed different expression patterns. Its expression was significantly increased in the breeding season and after E 2 treatments, but not as dramatically as erα. The ER antagonists, especially the Cyclofenil (a specific antagonist of ERβ), significantly attenuated the upregulation of erβ2 mRNA expression induced by E 2 . These data indicate that hepatic ERβ2 expression is regulated by estrogen. Interestingly, Cyclofenil also downregulated the increased expression of the erα and vtgs induced by E 2 , but DPN upregulated their expression. As the hepatic erβ1 expression was downregulated by E 2 and was not sensitive to Cyclofenil, the actions of Cyclofenil and DPN on the expression of erα and vtgs may possibly occur via ERβ2. To verify whether ERβ2 could regulate the expression of the erα and vtgs, we established a transfection system that transfers the siRNA into primary grouper hepatocytes using an electroporation method. Two siRNAs against erβ2 were employed to specifically knock down the hepatic erβ2 expression, and a significant decrease in the E 2 -induced expression of erα and vtgs in the hepatocytes was also observed. However, unlike ICI182780, neither Cyclofenil nor siRNA knockdown could fully abolish the E 2 induction of erα and vtgs expression, and the DPN-mediated upregulation of erα and vtgs expression was not as drastic as the E 2 treatment, suggesting that ERβ2 may partially contribute to the estrogenic regulation of erα and vtgs expression in grouper.
The results above showed that ERβ2 could promote the hepatic erα expression by mediating the estrogen action, but the molecular mechanism remains unclear. At present, it is clear that estrogen's regulation of downstream gene expression occurs through the E 2 -ER-ERE pathway. The E 2 -ER dimer complexes act on the ERE site in the promoter to initiate or repress gene expression [27,28]. Studies have shown that there are several ERE sites in the erα promoter region of rainbow trout and zebrafish [28]. However, there is no complete ERE site in grouper' erα promoter, but three 1/2 ERE sites, and some transcription binding sites that could mediate the estrogen signal, including one Sp1 binding site, nine Ap1 binding sites, and one CRBE binding site, were observed. Subsequent experiments demonstrated that E 2 significantly enhanced the grouper erα promoter activities in the presence of ERβ2 in HEK293T cells. Deletion and site-directed mutagenesis analysis indicated that the ERE#2 site is required to maintain the E 2 induction of erα promoter activities, suggesting that ERβ2 could modulate erα transcriptional activity through a half-ERE-dependent classical mechanism in grouper. Our previous study on grouper has shown that the half-ERE site in the gene promotor could initiate the estrogen response [29].
The present data suggest that ERβ2 is involved, directly and indirectly, in the regulation of VTG production in grouper. As promoters of vtg genes contain E 2 response sites, ERβ2 may bind to these sites and initiate the vtg expression in the hepatocytes directly. In addition, some studies proposed that ERα and ERβ2 may interact as heterodimers to drive the vtg expression [7]. However, this hypothesis requires further evidence in fish species. Indirectly, E 2 -activated ERβ2 enhances the hepatic erα expression, and then ERα promotes VTG synthesis, as our results demonstrate that E 2 could regulate the erα expression in the presence of ERβ2. During the breeding season, VTG is produced rapidly and efficiently in the female bloodstocks. The collaboration between ERα and ERβ subtypes may accelerate this process. In goldfish, a model for E 2 -mediated vitellogenesis was proposed: both ERβ subtypes can induce some vitellogenesis and are necessary for the normal E 2 induction of erα expression. The increased erα expression would sensitize the hepatocyte to further E 2 stimulation, and thus promote VTG production in the breeding season. This model may not be entirely suitable for grouper. Our study shows that grouper ERβ1 may not be involved in the regulation of VTG production. Although blocking ERβ2 action and the knockdown of its expression could result in the reduction of erα expression, whether ERβ2 is required for the E 2 -mediated upregulation of erα expression is uncertain and awaits further study. Therefore, the roles of ER subtypes, especially ERβ receptors, in mediating estrogen action may be different across fish species, suggesting that fish have evolved different mechanisms to regulate VTG synthesis.
In summary, we investigated the roles of ER subtypes towards VTG production in orange-spotted grouper. The results revealed that the physiological actions of estrogen on VTG production are primarily mediated by ERα, and that ERβ2 may play a role in promoting the E 2 -mediated upregulation of erα and vtg expression. This study will lead to a better understanding of vitellogenesis and paves the way for improving the oocyte quality of groupers.

Animal and Sample Collection
Eighteen adult orange-spotted groupers (5-years old; body weight: 3.21~4.46 kg; standard length: 46.93~60.16 cm) were obtained from the Hainan Chenhai Fisheries Co., Ltd., in Sanya (Hainan Province, China). The gonadal status of fish was extracted with a catheter and syringe, immediately fixed in Bouin's fluid, and checked by histological examination. The methods were carried out as described by Tang et al. [30]. Livers from non-breeding females (n = 6) were collected in February 2021, and livers from males (n = 6) and breeding females (n = 6) were collected in April 2021. Fish were anesthetized with MS-222 (Sigma) and dissected. The samples were frozen in liquid nitrogen and stored at −80 • C for RNA extraction. Juvenile female groupers (1.5-years old; body weight: 0.70~1.01 kg; standard length: 18.35~22.18 cm) used for intraperitoneal injection and primary hepatocyte culture were provided by Yuxiangzi Aquatic Industry Co., Ltd. in Yangjiang (Guangdong Province, China). All animal experiments were conducted in accordance with the guidelines and approval of the Animal Research and Ethics Committees of the Sun Yat-sen University.

Intraperitoneal Injection Experiment
Sixteen healthy juvenile female groupers (1.5-years old; body weight: 0.70~1.01 kg; standard length:18.35~22.18 cm) were randomly divided into two groups, and kept in two 3000 L tanks with continuous seawater and air supply. Groupers were acclimatized for a week before the experiment. The fish were exposed to natural photoperiods and water temperature (25.5-28.2 • C) and fed with commercial feed once a day. E 2 was dissolved in DMSO, and then diluted in physiological saline solution. The fish in experimental group were intraperitoneally injected with E 2 solution (5 mg/kg body weight). The physiological saline injected group was used as control. The livers were collected at 24 h after injection. The mRNA levels of vtgs and ers were quantified by Quantitative real-time PCR. The content of VTG protein in the liver was measured by using a VTG ELISA kit (Cusabio, Wuhan, China) following manufacturer's protocol.

Isolation, Primary Culture, and Static Incubation of Grouper Hepatocytes
Livers from six juvenile female groupers were excised and washed three times in ice-cold Ca 2+ /Mg 2+ -free HBSS (Gibco, Carlsbad, CA, USA) containing 1% antibioticantimycotic (Gibco). The tissue fragments were further diced to about 1 mm in thickness using a razor blade and washed three times in Ca 2+ /Mg 2+ -free HBSS containing 1% antibiotic-antimycotic at room temperature. Subsequently, slices were digested with trypsin-EDTA (0.25%) (Gibco) at 28 • C for 20 min, then were mechanically dispersed into individual cells by gentle pipetting and filtered through a 30 µm sterilized nylon mesh. Cells were harvested by centrifugation (30× g for 5 min, twice) and were resuspended in L-15 medium (Gibco) containing 1% antibiotic-antimycotic and 10% FBS (Gibco). The cell viability test and counting were performed using a trypan blue exclusion test. Only cells with >95% viability were used for subsequent experiments. The number of hepatocytes was adjusted to approximately 1.5 × 10 6 cells/mL/well and seeded as a monolayer into 24 well-plates precoated with 5 µg/mL PEI (Sigma, Saint Louis, MO, USA). After preincubation at 28 • C for 24 h, the medium was removed and replaced with a fresh medium containing the test drugs (E 2 , ER antagonists and agonist). In the E 2 treatment experiment, the doses of E 2 were 0.001, 0.01, 0.1, and 1 µM. In ER antagonists' experiments, the dose of E 2 was 0.01 µM, and 1 µM of ICI182780 or Cyclofenil was used. In ER agonist's experiment, 1 µM DPN was added into the medium. Each experiment had three groups (4 wells per group). After 24 h incubation, the culture media were collected for VTG protein assay, and cells were harvested for gene expression analysis.

Knockdown Experiment
Small interfering RNAs (siRNAs) ( Table 1) were synthesized by GenePharma according to the sequence of erβ2 (GenBank Accession No GU721078.1). Transfection experiment was performed on an electroporation instrument NEPA21 (NEPAGENE, Chiba, Japan). Primary hepatocytes (1.5 × 10 7 cells per mL) were suspended in the opti-MEM and mixed with siRNA at room temperature for 10 min. The mixture was transferred into electrode chamber and subjected to electroporation. The program was 195 V, 2.5 ms, 2 times, and 50 ms interval for poring pulse, and 20 V, 50 ms, 5 times, and a 50 ms interval for transfer pulse. After electroporation, the mixture was immediately transferred to L-15 medium containing 10% FBS in a 24-well plate. After preincubation at 28 • C for 24 h, the medium was removed and replaced with a fresh medium containing E 2 . After 24 h incubation, cells were harvested for gene expression analysis. The electroporation experiment was performed in triplicate. To analyze the transfection efficiency, siRNA-NC-Cy 3 (Cy 3 -labeled negative control siRNA) was transfected into primary hepatocytes by electroporation according to the protocol above. After 24 h incubation, the fluorescence value was detected by Beckman Cytoflex flow cytometry (Beckman Coulter Inc., Brea, CA, USA), and analyzed by CytExpert software (Beckman Coulter Inc., Brea, CA, USA).

Quantitative Real-Time PCR
To examine the mRNA levels of vtgs (vtgaa, vtgab, and vtgc) and ers (erα, erβ1, and erβ2), total RNA was extracted with Trizol reagent (Invitrogen, Carlsbad, CA, USA) and was reverse transcribed with a ReverTra Ace qPCR-RT Kit (Toyobo, Osaka, Japan) following the manufacturer's protocol. Quantitative real-time PCR reactions were performed on a Roche Light-Cycler 480 real-time PCR system using SYBR Green I Master (Roche, Basel, Switzerland) according to the manufacturer's protocol. Briefly, 500 ng template, 5 µL SYBR Green I Master, and 0.2 µL forward (reverse) primer (10 mM) were mixed in a tube and supplemented with water to 10 µL. The program was as follows: denaturation at 95 • C for 10 min, followed by 40 amplification cycles of 95 • C for 10 s, 58 • C for 20 s, and 72 • C for 20 s. After amplification, melting curve analysis was carried out to confirm the specificity of amplification. Fluorescence signals were converted to threshold cycle (Ct) values. Beta-actin was used as reference gene to normalize the expression values. The relative mRNA levels were calculated by using 2 −∆∆CT method. The sequences of gene-specific primers are listed in Table 2.

Gene
Primer Name 5 to 3 Sequence

Plasmid Construction and In Silico Analysis of Promoter
Using a SMARTer RACE cDNA Amplification Kit (Takara, Osaka, Japan), the transcriptional start site (TSS) of erα gene were confirmed. The 5 -flanking region (1434 bp) upstream of the transcriptional start site was cloned from the genome of orange-spotted grouper. The transcription factor binding sites in the promoter region were predicted using a virtual database online (http://alggen.lsi.upc.es/cgi-bin/promo_v3/promo/promoinit. cgi?dirDB=TF_8.3 (accessed on 2 July 2021)). Based on the sequence of 5 -flanking region, different erα promoter deletion fragments (1057 bp, 825 bp and 555 bp) upstream of TSS were cloned. These fragments were digested with endonuclease XhoI and KPNI (New England Biolabs, Ipswich, MA, USA) and then inserted into pGL4.10 vector (Promega, Madison, WI, USA). Based on the deletion analysis, three transcription factor binding sites of erα promoter were mutated using Mut Express II Fast Mutagenesis Kits (Vazyme, Nanjing, China) according to the manufacturer's protocol, and subsequently subcloned into pGL4.10 vector. In addition, the open reading frame (ORF) of erβ2 was cloned and subcloned into pcDNA3.1. All constructs were verified by sequencing to ensure there was no mismatch. Primers used for rapid amplification, site-directed mutagenesis, and plasmid construction are listed in Table 3.

Dual-Luciferase Assay
Human embryonic kidney cells, (HEK293T) purchased from Shanghai Institute of Biochemistry and Cell Biology (SIBCB), were seeded on a 48-well plate and cultured in phenol red-free DMEM medium (HyClone, Logan, UT, USA) containing 5% FBS and 1% penicillin and streptomycin at 37 • C with 5% CO 2 overnight. Cells in each well were co-transfected with 250 ng of promoter plasmid vector, 100 ng of ERβ2 expression vector, and 10 ng of pRL-TK vector (Promega). Six hours after transfection, the medium was replaced with fresh DMEM medium containing different doses of E 2 . After 24 h incubation, cells were harvested and lysed by cell lysis buffer (Vazyme). The luciferase activity was measured by the Dual-Luciferase Report Assay Kit (Promega).

Statistical Analysis
All statistical analyses were performed using SPSS 23.0. All data are presented as means ± SEM. Statistical differences were estimated by unpaired Student's t-test or oneway ANOVA followed by Tukey's post-hoc test. Values were considered significantly different at p < 0.05.