s2fieldsofstudy: Biology Medicine

text:
G-Protein-Coupled Estrogen Receptor Expression in Rat Uterine Artery Is Increased by Pregnancy and Induces Dilation in a Ca2+ and ERK1/2 Dependent Manner

Increasing levels of estrogens across gestation are partly responsible for the physiological adaptations of the maternal vasculature to pregnancy. The G protein-coupled estrogen receptor (GPER) mediates acute vasorelaxing effects in the uterine vasculature, which may contribute to the regulation of uteroplacental blood flow. The aim of this study was to investigate whether GPER expression and vasorelaxation may occur following pregnancy. Elucidation of the functional signalling involved was also investigated. Radial uterine and third-order mesenteric arteries were isolated from non-pregnant (NP) and pregnant rats (P). GPER mRNA levels were determined and—concentration–response curve to the GPER-specific agonist, G1 (10−10–10−6 M), was assessed in arteries pre-constricted with phenylephrine. In uterine arteries, GPER mRNA expression was significantly increased and vasorelaxation to G1 was significantly enhanced in P compared with NP rats. Meanwhile, in mesenteric arteries, there was a similar order of magnitude in NP and P rats. Inhibition of L-type calcium channels and extracellular signal-regulated kinases 1/2 significantly reduced vasorelaxation triggered by G1 in uterine arteries. Increased GPER expression and GPER-mediated vasorelaxation are associated with the advancement of gestation in uterine arteries. The modulation of GPER is exclusive to uterine arteries, thus suggesting a physiological contribution of GPER toward the regulation of uteroplacental blood flow during pregnancy.

Introduction
Normal pregnancy is associated with physiological adaptations of the maternal vasculature and involves significant changes to both systemic and tissue-specific vascular beds [1,2]. These adaptations are widely triggered by a decrease in vascular tone and an increase in uterine blood flow, which is essential to sustaining sufficient uteroplacental perfusion required by the developing fetus throughout pregnancy [3,4].

GPER mRNA Levels Change following Gestational Age in Uterine Arteries
In uterine arteries, mRNA levels of GPER were significantly higher in P7 rats, and this increase was even greater in P14 compared with NP rats (p < 0.05 P7 and P14 vs. NP; Figure 1A). A trend in the increase of GPER expression was assessed by qPCR in uterine arteries following pregnancy progression (p= 0.06, P7 vs. P14; Figure 1A). In contrast, in mesenteric arteries, there were no differences in the mRNA levels of GPER between pregnant and non-pregnant rats or throughout gestation ( Figure 1B). has been shown that natural estrogens stimulate the production of NO and other endo thelium-derived factors [13][14][15][16], and L-type Ca 2+ channels may be a target of acute estro genic vasodilatory effects during pregnancy [17]. Estrogen-induced vasodilation has the beneficial role of maintaining arterial health [18,19], with a potentially different degree of vascular effects according to tissue-specifi vascular beds. In particular, estrogens exhibit potent vasodilatory properties reaching their greatest effects on the reproductive tissues, where systemic infusion of these hor mones increases uterine blood flow up to 10-fold in non-pregnant ewes [20].
Previous studies have established that estrogens exert a functional role in the vascu lar adaptation to pregnancy and regulation of uterine blood flow [21,22] through binding to the estrogen receptor (ER), Erα, and ERβ [16,23,24]. In addition to the action elicited through these receptors, estrogens activate the G protein-coupled estrogen recepto (GPER) signalling leading to acute vascular effects [25,26]. Notably, we have demon strated that GPER is expressed in rat uterine arteries, and its activation by the specifi agonist G1 induces endothelium-dependent vasorelaxation through the NO-cyclic guano sine monophosphate (cGMP) pathway [27]. Moreover, previous evidence has shown tha -plasma levels of estrogens increase progressively in human pregnancy [28] as well as in pregnant rats [29].
However, it is currently unknown whether functional GPER expression changes dur ing pregnancy in the uterine and systemic circulation. Here, we ascertained that preg nancy modulates GPER expression and GPER-mediated vasorelaxation only in uterin arteries. In addition, we investigated mechanisms underlying GPER vasorelaxation.

GPER mRNA Levels Change following Gestational Age in Uterine Arteries
In uterine arteries, mRNA levels of GPER were significantly higher in P7 rats, and this increase was even greater in P14 compared with NP rats (p < 0.05 P7 and P14 vs. NP Figure 1A). A trend in the increase of GPER expression was assessed by qPCR in uterin arteries following pregnancy progression (p= 0.06, P7 vs. P14; Figure 1A). In contrast, in mesenteric arteries, there were no differences in the mRNA levels of GPER between preg nant and non-pregnant rats or throughout gestation ( Figure 1B). Pregnancy increases GPER expression in uterine but not in mesenteric arteries. GPER ex pression was quantified by qPCR in uterine (A) and mesenteric (B) arteries from non-pregnant (NP and pregnant rats at 7 (P7) and 14 (P14) days of gestation. GPER mRNA was significantly higher in uterine arteries from both P7 and P14 rats compared with NP rats (A). Expression levels were no affected by pregnancy in mesenteric arteries (B). * p < 0.05, P7, P14 vs. NP, by one-way ANOVA. = 3 independent experiments per group performed in triplicate.

Figure 1.
Pregnancy increases GPER expression in uterine but not in mesenteric arteries. GPER expression was quantified by qPCR in uterine (A) and mesenteric (B) arteries from non-pregnant (NP) and pregnant rats at 7 (P7) and 14 (P14) days of gestation. GPER mRNA was significantly higher in uterine arteries from both P7 and P14 rats compared with NP rats (A). Expression levels were not affected by pregnancy in mesenteric arteries (B). * p < 0.05, P7, P14 vs. NP, by one-way ANOVA. n = 3 independent experiments per group performed in triplicate.

Pregnancy Modulates Vascular Reactivity Response to G1 in Uterine but Not in Mesenteric Arteries
To provide further insights into the role of GPER, we then ascertained that the GPER agonist, G1, induces vasorelaxation of uterine ( Figure 2A) and mesenteric ( Figure 2B) arteries in a dose-dependent manner in both pregnant and non-pregnant rats. Notably, in uterine arteries we found a significant effect of pregnancy on G1-mediated vasorelaxation; it was modulated by gestational age (Figure 2A). G1-vasorelaxation was significantly higher in P7 compared with NP rats, with a maximum efficacy of 82.4 ± 6.0% in P7 vs. 66.5 ± 3.7% in NP rats ( Figure 2A). The effect of G1 on uterine arteries significantly increased following the progression of pregnancy, with a maximum vasorelaxation of 97.8 ± 2.5% in P14 rats ( Figure 2A) at 10 −6 M concentration of G1 ( Figure 2A). In contrast, in mesenteric arteries, dose-responses to G1 were similar between pregnant and non-pregnant rats and throughout gestation ( Figure 2B).

Pregnancy Modulates Vascular Reactivity Response to G1 in Uterine but Not in Mesenteric Arteries
To provide further insights into the role of GPER, we then ascertained that the GPER agonist, G1, induces vasorelaxation of uterine ( Figure 2A) and mesenteric ( Figure 2B) arteries in a dose-dependent manner in both pregnant and non-pregnant rats. Notably, in uterine arteries we found a significant effect of pregnancy on G1-mediated vasorelaxation; it was modulated by gestational age (Figure 2A). G1-vasorelaxation was significantly higher in P7 compared with NP rats, with a maximum efficacy of 82.4 ± 6.0% in P7 vs. 66.5 ± 3.7% in NP rats ( Figure 2A). The effect of G1 on uterine arteries significantly increased following the progression of pregnancy, with a maximum vasorelaxation of 97.8 ± 2.5% in P14 rats ( Figure 2A) at 10 −6 M concentration of G1 ( Figure 2A). In contrast, in mesenteric arteries, dose-responses to G1 were similar between pregnant and non-pregnant rats and throughout gestation ( Figure 2B).

Figure 2.
Pregnancy modulates concentration-responses to G1 in uterine arteries. In uterine arteries (A) of pregnant rats, vasorelaxation to G1 was significantly higher at 7 days (P7, n = 8) compared with non-pregnant (NP, n = 8) rats, and this response was even greater at 14 days of pregnancy (P14, n = 8).

Discussion
The present study demonstrates that (1) expression of GPER is increased, (2) vasorelaxation to G1 is enhanced as pregnancy progresses; GPER-mediated vasorelaxation oc-

Discussion
The present study demonstrates that (1) expression of GPER is increased, (2) vasorelaxation to G1 is enhanced as pregnancy progresses; GPER-mediated vasorelaxation occurs through (3) L-type Ca 2+ channels and (4) ERK1/2 pathway in uterine arteries of pregnant rats; (5) functional expression of GPER does not change through gestation in mesenteric arteries. Activation of GPER has been shown to cause vasorelaxation in rat aorta, carotid, mesenteric, renal, cerebral, and coronary arteries [30][31][32][33][34]. In addition, sexrelated controversial effects and stark differences between vascular beds have also been reported [32,[35][36][37][38]. Although considerable effort has been made to appreciate the role of GPER in the vasculature in both human and animal models (Supplementary Table S1, in Supplementary Materials), more studies are needed to fully delineate the physiological importance of GPER.
GPER together with the cognate receptors, ERα and ERβ, mediates physiological functions of estrogens, including the crucial role of regulating vascular function [25,26]. The agonist-receptor-coupling with the consequent production and release of NO and other endothelium-derived factors, determine the degree to which extent estrogens may elicit vasorelaxation [39,40]. In the present investigation, pregnancy significantly increased GPER-mediated vasorelaxation in uterine arteries, and this effect was modulated by gestational age.
To elucidate whether the modulating role of pregnancy on GPER-mediated vasorelaxation may be vascular-bed specific, we applied G1 in a cumulative fashion on arteries obtained from the mesenteric circulation. In agreement with other studies [30,34], G1 elicited a vasorelaxing effect in mesenteric arteries. Contrary to our findings in the uterine vasculature, the similarity of responses to G1 (41-52% maximum vasorelaxation) in mesenteric arteries of both -NP and P rats, as well as at different gestational periods, suggests that pregnancy-induced modulation of GPER vascular function is specific of the uterine circulation.
It may be implied that pregnancy-associated changes of the uterine reactivity to G1, could predict an up-regulation of GPER expression in this vascular bed. In fact, evaluation of GPER expression revealed that pregnancy modulates mRNA levels of GPER in uterine arteries, and a trend towards increase in GPER expression was observed as pregnancy progressed. A limitation of this study is that we investigated only mRNA levels of GPER. Further additional data in terms of protein expression are needed to fully support our conclusion. However, we may speculate that an increase in mRNA corresponds with an increase in the relative protein expression of GPER. Our results demonstrated a significant increase in GPER mRNA levels already after 7 days of pregnancy, with GPER expression exceeding a one-fold increase in P14 compared with NP rats. Contrariwise, we found no changes in GPER expression in mesenteric arteries, which could explain our functional results, and suggests again that pregnancy-induced modulation of GPER is specific to uterine arteries.
The lack of modulation in mesenteric arteries may be attributed to the different adaptation of the systemic circulation to pregnancy, in terms of both reactivity and remodelling [2], with the reproductive system undertaking a greater role in both maintaining vascular adaptation [5,24] and regulating uterine blood flow throughout gestation [22]. It is also plausible that the action of estrogens is more powerful at the uterine site than in the systemic circulation through the upregulation of GPER in uterine arteries, which in turn makes this vasculature more sensitive to estrogens. An overexpression of GPER, due to increased sensitivity to estrogens, has been observed in the human melanoma tissues [41].
Although we have previously demonstrated that GPER activation elicits endotheliumdependent vasorelaxation via the NO-cGMP pathway in rat uterine arteries [27], our further elucidation of the functional signalling identified a possible smooth muscle-related mechanism involved in the vascular responses. In our current study, inhibition of L-type Ca 2+ channels with verapamil caused about a three-times reduction of the vasorelaxation, thus suggesting that L-type Ca 2+ channels participate in the GPER-mediated vasorelaxing effect. This conclusion is supported by other studies that have provided evidence that GPER mediates rapid cell signalling via stimulation of intracellular Ca 2+ mobilisation [26,[42][43][44][45]. It has been demonstrated that extracellular application of G1 on human aortic vascular smooth muscle cells (VSMCs) increases intracellular Ca 2+ concentrations, and this effect is attenuated by GPER silencing [43]. Pre-incubation with the L-type Ca 2+ channel agonist Bay K8644, enhanced the vasodilatory estrogenic effects on placental arteries [17]. Furthermore, infusion of estradiol increased perfusion pressure and coronary resistance through a non-genomic molecular mechanism involving activation of the L-type Ca 2+ channels in rats [46]. Additionally, inhibition of L-type Ca 2+ channels with nifedipine significantly reduced G1-induced intracellular Ca 2+ increase in myometrial cells [47], and activation of GPER inhibited the endothelin-1-stimulated increase of intracellular Ca 2+ concentrations in VSMCs of murine carotid arteries [45]. Although there is heterogeneity amongst cells on the mechanism through which GPER induces Ca 2+ mobilisation [47], both the uptake of extracellular Ca 2+ and the release of the ion from intracellular Ca 2+ stores in cells may be involved in the GPER molecular signalling.
Studies conducted mainly in cell lines have shown that either estrogens or the GPERselective agonist G1, can activate a number of protein kinases, including ERK1/2, via a non-ER dependent mechanism that requires increased intracellular concentrations of Ca 2+ [48][49][50][51][52] through the L-type Ca 2+ channels [51]. It has been demonstrated that ERK 1/2 is a signalling component of the GPER cascade that induces NO production in human endothelial cells [53], with the same mechanism [54] through which estrogen-mediated rapid intracellular signalling regulates NO synthesis [55].
Therefore, once the participation of the L-type Ca 2+ channels was established, we sought to investigate the involvement of ERK 1/2 signalling in the vascular effects elicited by GPER in uterine arteries.
Incubation with the ERK 1/2 inhibitor, PD 098059, showed that vasorelaxation in response to G1 was attenuated by 24%, which is suggestive of a partial contribution of ERK1/2 in the mechanism of action of GPER in uterine arteries. The partial inhibition suggests that other protein kinases may be involved in the signalling cascade activated by GPER and warrant further investigation.
Consistent with our data, other studies provided evidence that GPER promotes activation of ERK1/2 in cancer cell lines [50,56], in rat heart [57,58], and in human umbilical vein smooth muscle cells [43].
In conclusion, we have shown that activation of GPER induces vasorelaxation of reproductive (uterine) and systemic (mesenteric) arteries. In uterine arteries, GPER-mediated vasorelaxation increases as pregnancy progresses. This can be explained by an upregulation of GPER gene expression in uterine arteries during pregnancy, which may be of importance in regulating hemodynamic changes of pregnancy in the reproductive system. Moreover, the emerging contribution of L-type Ca 2+ channels and ERK 1/2 signalling provides a new perspective for understanding the vascular mechanisms involved in GPER-mediated vasorelaxation of uterine arteries in pregnant rats. The modulating effect of both GPER expression and GPER vasorelaxation, associated with gestational age, is supportive of a physiological role for GPER in the uterine vascular adaptation and may offer a novel therapeutic target through which to selectively improve local uterine circulation in pregnancies with compromised uteroplacental blood flow.

Experimental Animals
The present study was performed in accordance with the European Guidelines on the protection of animals used for scientific purposes (

Pressure Myography
Radial uterine and third-order mesenteric arteries (diameter < 300 µm) obtained from NP, P7, and P14 rats were dissected free from perivascular connective and adipose tissue in ice-cold HEPES-PSS. Arterial segments (1-2 mm long) were transferred to the chamber of a small-vessel pressure myograph. One end of the artery was tied onto a glass cannula, any luminal content was flushed off by increasing intraluminal pressure before securing the distal end onto a second glass cannula using a servo-null pressure system (Living Systems Instrumentation, St. Albans City, VT, USA). All arteries were continuously superfused with HEPES-PSS at 37 • C, pressurised to 50 mmHg, and equilibrated for 45 min under no-flow conditions before the experiment started. Lumen diameter was measured by trans-illuminating each arterial segment in conjunction with data-acquisition software (Ionoptix, Westwood, MA, USA) to continuously record lumen diameter. Following equilibration, arteries were pre-constricted with phenylephrine (10 −7 ) to produce a 40-50% reduction in baseline diameter [59]. Once pre-constriction was achieved and remained stable for about 10 min, a concentration-response curve to the specific agonist of GPER, 1-(4-(-6-Bromobenzol(1,3)diodo5-yl)3a,4,5,9b-tetrahidro-3Hcyclopenta(c-)ethenone-8yl)ethenone (G1; 10 −10 -10 −6 M), was obtained (a representative pressure myography trace is reported in Supplementary Materials as Figure S1). To investigate the contribution of L-type calcium channels and ERK1/2 signalling, uterine arteries were pre-incubated with the inhibitors verapamil (10 −5 M) or PD 098059 (10 −6 M), respectively, for 30 min before pre-constriction with phenylephrine (10 −6 M) to produce a 40-50% reduction in baseline diameter, and vasorelaxation to G1 was assessed. HEPES-PSS without Ca 2+ plus papaverine (10 −4 M) was added to calculate the maximum vessel diameter at the end of each experiment.

Total RNA Extraction and Quantitative Real-Time Polymerase Chain Reaction (qPCR)
Uterine and mesenteric arteries were isolated as described above and preserved at −80 • C. Frozen samples were homogenised with a motor-driven homogeniser, and total RNA was isolated using Trizol reagent (Life Technologies, Milan, Italy), according to the manufacturer's instructions. RNA was quantified spectrophotometrically, and quality was checked by electrophoresis, through agarose gels stained with ethidium bromide. Total cDNA was synthesised from RNA by reverse transcription using the murine leukemia virus reverse transcriptase (Life Technologies, Milan, Italy), following the protocol provided by the manufacturer. The expression of GPER was quantified by real-time qPCR using the Step One TM sequence detection system (Applied Biosystems Inc., Milan, Italy). Specific primers for GAPDH and GPER were designed using Primer Express version 2.0 software (Applied Biosystems). The following sets of primers were used: for GAPDH (internal control), 5 -CAAGGCTGTGGGCAAGGT-3 (forward), 5 -GGAAGGCCATGCCAGTGA-3 (reverse); for GPER, 5 -CCTGGACGAGCAGTATTACGATAT-3 (forward), 5 -CGCTGCTGTACATGT TGATCTG-3 (reverse). Assays were performed in triplicate.

Drugs and Chemicals
G1 and PD098059 were purchased from TOCRIS (R&D Systems, Milan, Italy); Verapamil was purchased from Santa Cruz (Heidelberg, Germany). All the other chemicals were purchased from Sigma-Aldrich (Milan, Italy).
G1 stock solution was dissolved in DMSO, and small aliquots were frozen prior to use. All the other drugs tested were prepared daily. Phenylephrine and papaverine were dissolved in HEPES-PSS; verapamil and PD 098059 inhibitors were dissolved in DMSO. All drugs were kept on ice, in lightproof vials, and further diluted in HEPES-PSS as required.

Statistical Analysis
Vasorelaxation to G1 was expressed as percentage of maximally vasorelaxed diameter. Data are expressed as mean ± SEM, and n is the number of rats used. Differences between groups were determined by one-way ANOVA (for relative mRNA levels), and two-way ANOVA for repeated measures analysis, followed by Tukey's post hoc test where appropriate (for vasorelaxant responses). Differences in vasorelaxation following incubation with inhibitors were determined using paired Student's t-test, performed by Prism 8 software (GraphPad Software Inc., La Jolla, CA, USA). Differences were considered significant at p < 0.05.


Computational Analysis of the Potential Impact of MTC Complex Missenses SNPs Associated with Male Infertility

Meiotic chromosomes endure rapid prophase movements that ease the formation of interhomologue recombination intermediates that drive synapsis, crossing over, and segregation process. To generate these fast moves, the meiotic telomere complex (MTC) enables telomere-inner nuclear membrane attachment during meiotic prophase I and transfers cytoskeletal signals via another complex: the LINC complex. Furthermore, disruption or mutations of any of the MTC genes (TERB1, TERB2, and MAJIN) alters telomere association with the nuclear envelope leading to impairment of homologous pairing and synapsis, a meiotic arrest, and consequently to male infertility. To decipher the effect of TERB1, TERB2, and MAJIN missense mutations on protein structure, stability, and function, different bioinformatic tools were used in this study including VEP, Mutabind2, Haddock, Prodigy, Ligplot, ConSurf, DUET and MusiteDeep. In total, thirty mutations were predicted to be deleterious using VEP web server: seventeen for TERB1, eleven for TERB2, and two for MAJIN. All these single nucleotide polymorphisms were further analyzed and only 11 SNPs (W8R, G25R, P649A, I624T, C618R, F607V, S604G, C592Y, C592R, G187W, and R53C) were found to be the most damaging by at least six software tools and exert deleterious effect on the TERB1, TERB2, and MAJIN protein structures and likely functions. They revealed high conservation, less stability, and having a role in posttranslational modifications. This in silico approach provides information to gain further insights about variants that might affect stability, change binding affinity, and edit protein-protein interactions to facilitate their identification and functional characterization associated with male infertility.

Introduction
Telomeres have critical meiosis-specific functions in the reduction of chromosome numbers by assuring the integrity of the genome during meiosis. One such function is the attachment of telomeres to trans-nuclear envelope protein complexes that join telomeres to motor proteins in the cytoplasm. The active movement of telomeres and chromosomes during the first meiotic prophase is enabled by these trans-nuclear envelope linkages between telomeres and cytoplasmic motor proteins. Movements of chromosomes/telomeres facilitate the meiotic recombination process and ensure high fidelity pairing of homologous chromosomes which is a prerequisite for their correct segregation during the first meiotic division [1].
Meiosis is a gametogenesis-specific cell division that entails unique chromosomal regulations such as pairing, synapsis, and homologous chromosome recombination. These processes are ensured by the dynamic chromosomal rearrangements that happen during meiotic prophase I which is prolonged and more complicated than mitosis. Indeed, it can last many days and account for up to 90% of the total time spent in meiosis [2]. During meiotic prophase I, chromosomes undergo rapid movements that enhance the constitution of interhomologue recombination intermediates underlying synapsis, crossing over, and segregation processes [3] and ensure the pairing and eventual recombination of homologous chromosomes and the resolving of undesirable entanglements between nonhomologous partners [2].
The recently discovered tripartite complex comprising telomere repeats-binding bouquet formation proteins 1 and 2 (TERB1 and TERB2) and membrane-anchored junction protein (MAJIN) has been demonstrated to be involved in these processes by tethering telomere ends to the nuclear envelope (NE), transmitting cytoskeletal forces via the LINC complex to drive these rapid movements, and assembling of meiotic mammalian telomeres on the NE during prophase I [4].
In mammals, the meiotic telomere complex (MTC) reuses and integrates the functions of two complexes, Shelterin and LINC, which otherwise play distinct significant roles outside of meiosis [3]. The MTC complex is initially sequestered onto the NE by the membrane-binding activity of MAJIN and then assembles on the telomeres when they approach the nuclear periphery at meiotic entry. The assembly of TERB1, TERB2, and MAJIN onto telomeres is needed for the attachment of the LINC complex proteins (SUN1 and KASH5), and it is required for a stable telomere-NE attachment before the telomere movement ( Figure 1) [4].

Identification of Mammalian Shelterin
Complex. Mammalian telomeres contain a specific protein complex that accomplishes two functions: protecting the natural chromo-some end from all factors of the DNA damage response within the cell and negative regulation of telomerase by segregation of its telomeric DNA substrate. Both functions of this complex refer to the name shelterin [5]. The shelterin complex includes six telomere-specific proteins TRF1, TRF2 (telomeric repeat binding factors family proteins 1 and 2), POT1 (protection of telomeres 1), RAP1 (Repressor/activator protein1), TIN2 (TRF1-interacting protein 2), and TPP1 (TINT1/PTOP/PIP) that protect telomeres from degradation, inhibits unnecessary repair mechanisms, controls telomerase production, and is involved in cellular senescence and age-related pathologies [6,7].
Because of their TRF homology domains, the shelterin complex's subunits TRF1 and TRF2 perform a fundamental role by recognizing telomeric DNA, and by recruiting the other shelterin constituents onto telomeric DNA [8,9], where the reformed shelterin core of TRF2-TIN2-TPP1-POT1 has a stoichiometry of 2 : 1 : 1 : 1 in vitro [10]. During prophase, the MTC complex matures, allowing the shelterin/telosome complex to be released from telomeric DNA [7]. Shelterin acts as a telomere recognition protein for the meiotic telomere complex (MTC), (Figure 1) externally connecting it to the LINC complex and allowing them to work together to achieve cytoskeletal attachments of meiotic telomere ends around the nuclear envelope [3].  BioMed Research International containing 1) and KASH5 (Klarsicht/ANC-1/Syne/homology 5) proteins assemble into the mammalian meiosis specific LINC complex (linker of nucleoskeleton and cytoskeleton), providing the binding sites for telomeres on the inner surface of the nuclear envelope (NE) [11,12] and procuring the physical linkage to transfer cytoskeletal force transduction from the cytoplasm to the nucleus [13][14][15]. KASH proteins interact with the cytoskeleton across the outer nuclear membrane, whereas SUN proteins connect with nuclear lamin A and emerin across the inner nuclear membrane [16][17][18][19]. In meiotic prophase I, the nuclear lamina endures an important reorganization where LINC complexes link with the meiotic telomere complex (MTC), thanks to movements assured by SUN1 and KASH5 proteins, which interact with microtubules via dyneindynactin ( Figure 1) [16,20].

The Meiotic
Telomere Complex (MTC). Another protein complex consisting of TERB1 (telomere repeat binding bouquet formation protein 1), TERB2 (telomere repeat binding bouquet formation protein 2), and MAJIN (membrane anchored junction protein) has been demonstrated to build up a second actual linkage for telomere connection to the nuclear envelope by assuring meiotic telomere attachments in mammals [21,22]. The meiosis-specific telomere regulator TERB1 is a molecular scaffold that simultaneously interacts with SUN1 and meiotic cohesin subunit SA3 through its N-terminal ARM (armadillo) repeat domain and C-terminal Myb domain, respectively, providing telomere attachment to the inner nuclear membrane (INM) and driving the chromosome movement essential for homologous pairing and recombination [21,22].
Through a region surrounding its TERB2-binding site, TERB1 interacts directly with shelterin constituent TRF1, including a peptide interaction that resembles TIN2 binding to the TRF1 dimeric cleft [4,23]. Consequently, Long et al. presumed that the disruption of the TERB1-TRF1 interaction impairs telomeric localization of TERB1 and SUN1 in spermatocytes and that the TERB1-TRF1 interface is specific and required for both in vitro and in vivo linking of TERB1 to TRF1 [4].
Therefore, an interaction of the meiotic telomere complex with TRF1 seems to be transient in the cell. MAJIN-TERB2-TERB1 assembles on the nuclear envelope and endures TRF1-dependent recruitment of telomere ends during leptotene and zygotene [4,22]; CDK activity at that point triggers the relocation of TRF1 to flanking areas in pachytene, with MAJIN-TERB2-TERB1 remaining related with telomere ends [22].
Thus, in mice, knocking out the Terb1 gene disturbs the whole interaction network and affects homolog pairing, synapsis, and recombination, resulting in early spermatogenesis and oogenesis [21]. However, the relevance of each of the TERB1-mediated contacts and their molecular mechanisms in meiosis remain unknowns [4].
MAJIN is a putative transmembrane (TM) protein, located at the inner surface of the NE; it serves as an anchoring component for the inner nuclear membrane through a transmembrane helix at its C-terminus [22]. MAJIN has DNA-binding activity, suggesting that it is involved in the stability of telomere attachment to the inner membrane of the nucleus. MAJIN and TERB1 are physically linked by TERB2, which binds MAJIN through its C-terminus and TERB1 through its N-terminus [4,22]. Individual disturbance of meiotic telomere complex components causes meiotic arrest in mice, with loss of telomere attachments and chromosomal movements, failure of DNA double-strand break repair, and disrupted synapsis [21,22]. Accordingly, the disturbance of telomere-NE connection impairs the homologous searching that causes subsequent meiosis defects, including abolished telomere-NE attachment, aberrant homologous pairing, disordered synapsis, and infertility in both sexes [24]. Collectively, these results reveal a linear interaction network within the TERB1, TERB2, and MAJIN proteins that shape together a stable complex that performs a major role in regulating the recruitment of telomeres to the NE ( Figure 1) [21][22][23][24].
Lately, computational analysis is designed to more accurately predict the impacts of mutations on proteins interactions. By using modern bioinformatics tools as a systematic in silico approach, our study is aimed at identifying the deleterious SNPs in the TERB1, TERB2, and MAJIN genes and predicting their significant pathogenic impact on the functions, stability, and structures of TERB1-TERB2 complex and TERB2-MAJIN complex.

Materials and Methods
In this study, we carried out a systematic in silico approach using modern bioinformatics tools to predict deleterious SNPs in the TERB1, TERB2, and MAJIN genes and identify their significant pathogenic impact on the functions and structures of the TERB1, TERB2, and MAJIN proteins ( Figure 2).

Prediction of Mutations Effects.
To predict the impact of the retrieved SNPs, the Ensembl Variant Effect Predictor (VEP) was utilized. This tool provides an indication of the effect of the amino acid change using protein biophysical properties [25]. These data can help with the interpretation of protein variants with no associated phenotype or disease data by estimating how deleterious a given mutation may be on the functional status of resulting protein [25]. Scores and predictions are calculated for all possible amino acid substitutions using different algorithms as follows.
The Sorting Intolerant from Tolerant (SIFT) algorithm is predicting whether an amino acid substitution affects a protein function. It takes into account the position of the mutation as well as the type of amino acid change. SIFT assesses the probability that an amino acid at a particular location will be tolerated based on the most frequently tolerated 3 BioMed Research International amino acid. The substitution is predicted to be deleterious, if the normalized value is smaller than a cutoff. Scores range from 0 to 1. The smaller the score, the more likely the SNP has damaging effect [26][27][28].
Polymorphism phenotyping (PolyPhen2) is a tool that predicts the possible impact of amino acid substitutions on the human protein structure and function using structural and comparative evolutionary considerations. PolyPhen2 results are available for human proteins. It classifies the substitution as probably damaging (score = 1 : 0) and possibly damaging or benign (score = 0 : 0) [29].
Other pathogenicity predictor scores such as Mutation-Taster, MetaSVM, and MetaLR are available for human data via VEP plugins: MutationTaster is used for a quick assessment of the disease-causing potential of DNA sequence changes. Evolutionary conservation, splice-site changes, protein feature loss, and changes that could affect the amount of mRNA are all examined using MutationTaster, which predicts an alteration as one of the four possible types: "A" (disease causing automatic), "D" (disease-causing), "N" (polymorphism), or "P" (polymorphism automatic) [30,31]; MetaSVM (SVM (Support Vector Machine)) is a linear model that can be used to solve classification and regression problems. The theory behind SVM is simple: it finds and creates a line (hyperplane) that separates the data into different classes. Prediction using SVM approach in Ensembl database reveals "T" as tolerated or "D" as damaging or deleterious. The score cutoff between "D" and "T" is 0, and higher scores are more deleterious [32]; In MetaLR, to predict the deleteriousness of missense variants, logistic regression (LR) is used to combine nine independent variant deleteriousness scores with allele frequency information. Variants are categorized as "tolerated" or "damaging"; a score between 0 and 1 is also provided where variants with higher scores are more likely to be deleterious [32].

Protein
Sequence. The three-dimensional structure of each complex was downloaded from RCSB Protein Data-Bank (https://www.rcsb.org/): The PDB ID for TERB1-TERB2 complex is 6J07 and for TERB2-MAJIN complex is 6GNY ( Figure 3).

Prediction of SNP Effects on Protein-Protein Interactions.
In order to estimate the impacts of single mutations on protein-protein interactions (PPI), we performed computational analysis on the human crystal structure of TERB1-TERB2 and TERB2-MAJIN complexes using MutaBind2 [33]. This tool compares binding affinity after mutations to predict whether they stabilize or destabilize the PPI by determining the overall change in binding free energies (ΔΔG) by providing the following results: The ΔΔG bind (kcal mol -1 ) which is the predicted change in binding affinity induced by a mutation; in deleterious (yes or no), MutaBind2 server classifies a mutation as deleterious if ΔΔG ≥ 1:5 or ≤−1:5 kcal mol -1 , ΔΔE vdw is the change of van der Waals interaction energy upon a mutation, ΔΔG solv is the change of polar solvation energy upon mutation, and ΔΔG fold is the change of stability of protein complex upon mutation [33].

Mutagenesis and Energy
Minimization. The chains of each complex (A and B of the 6J07 complex and C and D of the 6GNY complex) were separated using the UCSF CHI-MERA program [34]. Mutant structures of TERB1, TERB2, and MAJIN proteins were created manually by PyMol software (ver.2.4, Schreödinger) [35]. Then, the energy minimization of structures was performed by YASARA Energy Minimization web server [36].  5 BioMed Research International 2.6. Molecular Docking Study. The High Ambiguity Driven protein-protein Docking (HADDOCK) web server was employed to perform protein docking. The HADDOCK2.4 uses chemical shift perturbation data resulting from NMR titration experiments, mutagenesis data, and bioinformatics predictions [37]. In the HADDOCK2.4 submission form, we docked the structures of the energy minimized mutated proteins with each other's and with the wild type; the protein and peptide structure files were uploaded as Molecules 1 and 2, respectively. The active residues that form the interaction of proteins with each other in each complex have been provided according to preliminary information from PDBsum which is one of several online databases that provide information on all experimentally determined structural models published by the Protein Data Bank (PDB) [38]. Rigid molecular docking with flexibility on both active protein side chains and peptide structures was executed using the default settings [38]. The docking score (haddock score) of each variant of TERB1, TERB2, and MAJIN is hereafter designated as the stability index of TERB1-TERB2 and TERB2-MAJIN complexes. Cluster with the lowest haddock score is the most likely conformation [39].

Binding Affinity Prediction.
Subsequently, the protein binding energy prediction (PRODIGY) web server was used to predict the binding affinity of the protein-protein complexes from their 3D structure [39,40]. Complexes with the lowest energy value (ΔΔG) have greater binding affinity and correlates with mutations that stabilize the protein structures [40].

Analysis of Protein-Protein
Interactions. Ligplot software uses the 3D coordinates of a protein and its bound protein/ ligand to automatically produce schematic diagrams. These diagrams show the pattern of interactions between the two molecules and are very useful for comparing different structures since they provide clear and helpful information of the intermolecular interactions and their strengths, including hydrogen bonds, hydrophobic interactions, and atom accessibilities [41].
2.9. Assessment of Protein Stability. The stability of the protein was checked using the Protein Stability Change Upon Mutation (DUET), which is an online server for an integrated computational approach to study missense mutations in proteins. It combines two complementary approaches (mCSM and SDM) in a consensus prediction, obtained by consolidating the results of the separate methods in an optimized predictor using Support Vector Machines (SVM) [42]. DUET predicts the change in stability by calculating  BioMed Research International the change in unfolding free energy (ΔΔG). It further defines whether these changes increase or decrease the stability of the protein. Positive ΔΔG value denotes that protein stability increased, and negative ΔΔG means that protein stability decreased [43]. Site Directed Mutator (SDM) is an algorithm of statistical potential energy function that calculates a stability score based on environment-specific amino acid substitution frequencies within homologous protein families, which is comparable to the free energy difference between wild-type and mutant proteins [44].
Mutation Cutoff Scanning Matrix (mCSM) relies on graph-based signatures used as a novel approach to the study of missense mutations and the prediction of their effects. These signatures encode distance patterns between atoms and are used to represent the protein residue environment and to train predictive models [45].

Phylogenetic Conservational Analysis.
To predict the evolutionary conservation of the amino acids in a protein sequence, we used a ConSurf bioinformatic tool that provides evolutionary profiles of each of the amino acids in the protein, based on phylogenetic relations between homologous sequences to reveal regions that are important for structure and/or function [46]. The tool also predicts the conservation score for each amino acid residue ranging from 1 to 9, where the score denotes the degree to which the amino acids are evolutionary conserved: 1-3 designate variable residues, 4-6 designate medium conserved scores, and 7-9 depict highly conserved residue [47,48].

Root Mean Square Deviation Calculation. YASARA
View is an open-source program for the molecular graphics, modeling, visualization, and analysis of the threedimensional protein's structures. The structural deviations 7 BioMed Research International between the native and mutated models were analyzed using the YASARA View program, by measuring the root mean square deviation (RMSD) which is the average distance between the atoms of the superimposed proteins [49]. Its values are considered to be reliable indicators of variability, where values superior than 0.15 were estimated to be significant and can affect protein function and/or structure [49].

Prediction of Posttranslational Modification Sites.
A posttranslational modification is crucial for cell signaling and affects the function of the protein. TERB1, TERB2, and MAJIN proteins' posttranslational modifications were predicted using MusiteDeep, an online service that includes a deep learning framework for predicting protein posttranslational modification (PTM) sites [50].

Retrieval of Deleterious SNPs.
Thirty missense SNPs were predicted to be deleterious using SIFT (where tolerance index score was ranged from 0 to 0.02). The SNPs predicted by SIFT were validated by PolyPhen2, MutationTaster, MetaSVM, and MetaLR. In particular, PolyPhen2 results showed that eight AA substitutions (R620C, I624T, S608C, F607V, R605Q, C592Y, C592R, and L189P) were predicted to be possibly damaging, while twenty-two SNPs were of PRODIGY showed that the lowest binding energy indicates greater binding affinity, while higher binding energy is giving less binding affinity ( Figure 5).
Based on heatmap results, regions with the lightest orange shade represent complexes with higher free binding energy (ΔΔG) which correlates with lower binding affinity.  (Table 3).
Mutated complex represented some differences in the hydrogen bonds and hydrophobic interactions: some contacts were added and others were removed. R620C-W8R complex showed that hydrogen bonds formed by Gly, Val, Ser (3), Arg (8), Leu, His, Leu, Ile, Thr, Asp, Asn, and Gln
S604G-G25V complex showed that hydrogen bonds constituted by Gly, Asp, Asn, Arg, Tyr, Gln, Thr, and His residues were added and Leu, Val, Asp, Glu, Thr, Asn, Arg, Gln, His, Lys, Ser, and Ile were removed. Some hydrophobic interactions involving Gln, Lys, Ser, Arg, Val, and Pro have been added, and Arg, Cys, Asn, Val, Lys, Glu, Pro, and His were eliminated.
The mutated complex C618R-F64S revealed that hydrogen bonds formed by Asp, Ile, Asn, Thr, Gly, Gln, Arg, and His amino acids were added, while those formed by Gly, Val, Ser, Arg, Lys, His, Leu, Ile, Ala, Thr, Glu, Asp, and Gln were canceled. For hydrophobic contacts, Lys, Gln, Pro, Val, and Arg amino acids were added and Arg, Cys, Asp, His, Pro, Asn, and Gly amino acids were deleted.
The mutated complex Y631N-G84D revealed that hydrogen bonds formed by Asp, Tyr, Asn, Arg, Cys, Ile, His, and Phe amino acids were added and hydrogen bonds formed by Leu, Val, Asp, Glu, Thr, Gly, Ser, Lys, Arg, and Ile amino acids were removed. For hydrophobic contacts,  (Table 3).
WT TERB1-G25R complex showed that hydrogen bonds formed by Asp, Gln, Tyr, Asn, Thr, Arg, Tyr, and His amino acids were added, while those formed by Val, Ser, Arg, His, Leu, Thr, Asp, Gln, Asn, and Ile amino acids were deleted. Five hydrophobic contacts (Ser, Asp, Glu, Arg, and Val) were added, and five (Arg, Asn, Glu, Pro, and Tyr) were removed.
The second wild complex TERB2-MAJIN undergoes 20 hydrogen bonds formed between amino acids (Ser, Lys (5) The mutated complex WT TERB2-R53H showed that six amino acids (Ser, Asp, Cys, Gly, Val, and Lys) were added to hydrogen bond interactions, while three (Arg, Gln, and Pro) were removed. The hydrophobic interactions listed two amino acids (Ile and Asp) as deleted and only one amino acid (Pro) as added.
Five amino acids (Arg, Ala, Gly, Ile, and Val) were added to hydrogen bonding of the mutated complex MAJIN-L189P, and five were removed (Pro, Gln, Ala, Ser, and Ile). Hydrophobic contacts added three new residues (His, Gln, and Pro) to the list, while four residues Leu, Ile, Val, and Gly were taken off ( Table 3). The mutated complex MAJIN-D191N showed that hydrogen bonds formed by Ala, Val, Asp, and Lys amino acids were dismissed, while two others were added (Leu and Asp). Some hydrophobic interactions implying Gly, His, and Asn have been added, and Asp, Val, and Ile have been eliminated (Table 3). About the hydrophobic contacts of the mutant complex L189P-R53C, two amino acids were added (Gln and Pro), and six were canceled (Ile, Asp, Gly, Glu, Leu, and Ser). In the hydrogen bonds, six amino acids (Asp, Gly, Cys, Glu, Gly, and Lys) were added, while one (Gln) was removed.
DUET, ConSurf, and MusiteDeep web servers were used to identify, respectively, the mutation effects on protein stability and the conserved domains in proteins and predict directly the posttranslational modification (PTM) site from the raw protein sequence. The mutated structures were also analyzed by YASARA View, with regard to evaluate conformational variations by calculating the RMSD.  (Table 4). Furthermore, other mutations increased the stability of the proteins by one or two algorithms (Table 4).

Phylogenetic Conservation Analysis of High-Risk SNP.
Compared to those in nonconserved regions, amino acids located in conserved regions were predicted to be highly damaging. ConSurf predicts amino acids to play structural or functional roles based on conservation and solvent accessibility. Residues are predicted as functional when they are 13 BioMed Research International highly conserved and exposed and as structural when they are highly conserved and buried.
Phylogenetic conservation analysis of TERB2 protein showed that G25, P90, G187, D191, and G198 are highly conserved residues and predicted to be exposed and functional. On the other hand, W8, F64, G84, and L189 are highly conserved residues and predicted to be buried and structural. Finally, R53 is highly conserved and predicted to be exposed and a functional residue in MAJIN protein (Figure 6(a)).
Analysis of TERB1 protein from the ConSurf server revealed that P649 and S604 are highly conserved residues and predicted to be exposed and functional, while I624, H621, C618, F607, and C592 are highly conserved residues and predicted to be buried and structural. K632, R620, S608, and R605 are exposed residues: R620 is variable, while K632 and S608 are conserved. Y631 is a conserved and buried residue (Figure 6(b)).

Root Mean Square Deviation (RMSD) Prediction.
In our study, 29 mutations out of 30 showed higher RMSD values which indicate greater variation between wild and mutant protein structures. Therefore, four variants C592R, G25V, G198R, and R53H reflected maximal structural dissimilarity with TERB1-TERB2 complex and TERB2-MAJIN complex (1.823, 1.4844, 1.4723, and 1.6674, respectively) and were therefore used for superimposing (Figure 7). One residue change G25R in TERB2 protein showed a RMSD value equal to 0 suggesting a minimal structural deviation between the native model and this mutation (Table 5).

Prediction of Posttranslational Modification Sites.
Musi-teDeep is applied to predict the effect of SNPs on posttranslational modification (PTM) process of the human TERB1, TERB2, and MAJIN proteins. MusiteDeep identified sites for O-linked glycosylation: R605 in TERB1 protein and W8, G84, G187, L189, D191, and G198 in TERB2 protein sequence (Table 6). S-Palmitoylation site was identified in C592 of TERB1 protein. No PTM site has been located on MAJIN protein sequence.
3.9. Cumulative Score (CS) Calculation. A cumulative score (CS) was calculated for the 30 damaging variants based on seven in silico tools used for predicting the effect of SNPs (SIFT, PolyPhen2, Mutabind2, Prodigy, stability by SDM, mCSM, DUET, and RMSD estimated by YASARA, and ConSurf) to further understand and justify how the most deleterious mutations were identified in this study ( Table 7).

Discussion
The results of several recently published papers have contributed to our understanding of the meiotic telomere complex assembly in mammals and its possible interaction with  showed that TERB1 = CCDC79 structure contains two armadillo repeats at the N terminus, a coiled coil motif in the middle region and, remarkably, a Myb domain at the C terminus, a structure which resembles that of TRF1. At this point, many roles were dedicated to TERB1 [21]. It is important for homologue pairing/synapsis and recombination and for male and female fertility. TERB1 assembles a large protein complex linking telomeres to the SUN1-KASH5 complex along the nuclear envelope in meiosis. Moreover, TERB1 binds to telomeres by forming a heterocomplex with TRF1, and its Myb domain performs a crucial role in enriching cohesin, thereby promoting the structural integrity of telomeres during prophase I [21]. In addition to TERB1, two downstream meiotic proteins, TERB2 and MAJIN, are recruited to the meiotic telomeres by TERB1 in mammalian meiosis [22]. Shibuya et al.
affirmed that MAJIN was the key player, named also after a Japanese word representing the "Genie in Aladdin's Lamp." MAJIN behaves as a membrane protein and sequesters telomere adaptors TERB1/2 to the INM. Thus, they concluded that TERB2-MAJIN plays a critical role in assembling SUN1-KASH5 at telomeres and achieving meiotic chromosome movement [22]. In 2017, Long et al. focused on the various key domains of TERB1 that mediate different protein interactions that are required for constructing meiosisspecific telomere structures that enable telomere attachment and movement along the NE for faithful homolog pairing in mammalian meiosis [4]. One year later, Dunce et al. reported the crystal structure of the MAJIN-TERB2 complex, revealing a 2 : 2 heterotetramer in which two TERB2 chains wrap around a core MAJIN globular dimer. This structure endures direct interaction with DNA and scaffolds assembly of the full 2 : 2 : 2 TERB1-TERB2-MAJIN of the meiotic telomere complex, which can recruit two TRF1 dimers. Together, this data confirmed that TRF1-mediated loading of telomeric DNA to the meiotic DNA leads to a
In 2019, Wang et al. concluded that the intact MTC complex is required for gametogenesis and fertility, and targeted disruption of the complex induces disordered synapsis and impairs meiotic double-strand break (DSB) repair and abolishes telomere-NE linkage, together leading to a complete meiotic arrest in prophase I in both sexes. Hence, both the MTC complex and SUN1 (the LINC complex) contribute to the stable telomere-NE association and are important for efficient progression of prophase I [24].
Last year, Salas-Huetos et al. strongly confirmed that disruption of any of the three genes coding for the meiotic telomere complex (MTC) can cause NOA in men. It has been demonstrated previously that disruption of Terb1 in the mouse results in meiotic arrest and impairment of homologous pairing and synapsis, ultimately resulting in infertility in both males and females [21]. Likewise, mice disrupted for either MAJIN or TERB2 display impaired synapsis, zygotene arrest, a lack of post-meiotic cells, and infertility phenotype as revealed by Salas-Huetos [51].
Taken together, they showed that the disruption or mutations of any of the MTC genes impairs telomere association with the inner nuclear membrane, which is requisite for chromosomal synapsis and chromosome movement during meiosis. Failure of these events triggers meiotic checkpoints leading to a meiotic arrest and consequently to male infertility/NOA [51]. The main goal of this work was exploring the effects of missense mutations on meiotic telomere complex proteins (TERB1, TERB2, and MAJIN) related to male infertility.
Bioinformatic analysis provides powerful insights into understanding and predicting the effects of mutations on protein structures, including how mutations alter protein stability, binding energy, binding affinity, evolutionary conservation and interactions with other proteins and their relation to susceptibility of genetic diseases. In the current study, various algorithms, namely, VEP tools (SIFT, PolyPhen2, and MutationTaster), Mutabind2, Haddock, Prodigy, Ligplot, DUET, ConSurf, and MusiteDeep, were used to identify the most deleterious SNPs of the TERB1, TERB2, and MAJIN genes. Accordingly, eleven SNPs out of thirty (W8R, G25R, P649A, I624T, C618R, F607V, S604G, C592Y, C592R, G187W, and R53C) scored the highest cumulative scores (CS = 6) predicted by different bioinformatic tools and were located in conserved regions with a score conservation of 9, which may affect the structure and function of proteins (Table 7).
Protein-protein docking aims to anticipate molecular complex's structure given the individual solved structures of its constituents, which can be determined experimentally or predicted. In the present study, docking was performed by HADDOCK web server that displays results as cluster's structures, where the first top cluster produced is deemed reliable and should be downloaded in order to be loaded on the PRODIGY web server to calculate binding energy. Binding affinity refers to the strength or interaction energy of the binding between the protein and its partner. According to the results of heatmap, the lighter colored areas repre-sent the complexes with higher binding energy and less binding affinity: R620C-W8R, S604G-W8R, S604G-G25V,  C618R-F64S, C618R-G84D, Y631N-G84D, and WT  TERB1-G25R of TERB1-TERB2 complex and WT TERB2-R53C, WT MAJIN-L189P, WT MAJIN-D191N, and L198P-R53C) of TERB2-MAJIN complex ( Figure 5). Generally, the steady state of a protein is the state with the lowest energy level. The energy score influences protein's stability and binding affinity which means: lower energy value correlates with a favorable mutation that stabilizes the protein structure and can potentially improve binding affinity, and higher energy values refer to mutations that decrease protein stability and potentially alter the binding affinity.
Missense mutations can affect protein-protein interactions via various molecular mechanisms, including changing folding free energy of interacting partners, modifying protein stability, or disrupting noncovalent interactions essential for complex formation and function. Binding affinity is also influenced by noncovalent interactions between biomolecules, such as the loss or gain of hydrogen bonding and/or hydrophobic contacts in the mutated complex compared to the wild-type complex. After visualization by Ligplot, we concluded that the eleven previous complexes with the lower binding affinity showed differences in protein-protein interactions where some contacts were added and others were deleted, which explains the diminution of binding affinity within these complexes (Table 3).
DUET online server was used to check the significant effect of mutations on the stability of protein structure. Out of the most deleterious variants, ten variants revealed negative unfolding free energy (ΔΔG) predicted by the 3 algorithms which indicate that they destabilized the protein's stability, except for the S604G that was destabilizing only by two algorithms ( Table 4).
The RMSD computation was used to assess the superposition when comparing native and mutated proteins. Among the eleven mutations, ten showed higher RMSD values, except for the residue change G25R that had a RMSD value = 0 indicating (Table 5).
Conservation analysis is important in unraveling whether the SNPs are found in a conserved region or not. Binding regions are known to be evolutionarily conserved, which has been evaluated in various studies to identify potential protein interaction interfaces. Thus, according to Miller and Kumar, highly conserved amino acids are located in biologically active sites. When these residues are substituted, biological activities are affected [52]. The results of the ConSurf analysis showed that the eleven amino acid changes (W8R, G25R, P649A, I624T, C618R, F607V, S604G, C592Y, C592R, G187W, and R53C) were located in conserved regions ( Figure 6).
MusiteDeep identified sites for O-linked glycosylation W8 and G187 in TERB2 protein sequence with a score of 0.864 and 0.771, respectively. Consequently, mutations at these positions might affect function and/or structure of the TERB2 protein (Table 6).
Most improvements to computational predictions of mutation effects have been achieved by the identification of SNPs and their impact on human health, which has been 16 BioMed Research International the most explored field in human genetics recently. In this study, we applied several tools to prioritize damaging SNPs and estimated the most deleterious ones based on binding energy and affinity evaluation, stability assessment, evolutionary conservation analysis, and posttranslational modification site prediction. However, in silico prediction of protein-protein interaction network in native and wild protein should be confirmed with extensive experiments and lab approaches to figure out the mechanism and impact of these mutations in susceptibility to male infertility.

Conclusion
In the current study, out of 32 missense mutations, 11 were cumulatively predicted to be high risk pathogenic SNPs (W8R, G25R, P649A, I624T, C618R, F607V, S604G, C592Y, C592R, G187W, and R53C) as they were agreed by six software. Combinations of multiple in silico tools provided information required to predict the effects of these mutations on the functional and structural deviations of TERB1, TERB2, and MAJIN proteins. This work has highlighted the need to take into account the type of proteinprotein interactions when characterizing the variants within the complexes. Hence, our result added to previous knowledge and supported earlier findings that may be helpful for further understanding the role of meiotic telomere complex (MTC) in male infertility. Nevertheless, bioinformatics tools cannot replace conclusive experiments, and their results should be verified by supplementary data.

Data Availability
All data are available from the corresponding author upon request.


Plastome characteristics and species identification of Chinese medicinal wintergreens (Gaultheria, Ericaceae)

Wintergreen oil is a folk medicine widely used in foods, pesticides, cosmetics and drugs. In China, nine out of 47 species within Gaultheria (Ericaceae) are traditionally used as Chinese medicinal wintergreens; however, phylogenetic approaches currently used to discriminating these species remain unsatisfactory. In this study, we sequenced and characterized plastomes from nine Chinese wintergreen species and identified candidate DNA barcoding regions for Gaultheria. Each Gaultheria plastome contained 110 unique genes (76 protein-coding, 30 tRNA, and four rRNA genes). Duplication of trnfM, rps14, and rpl23 genes were detected, while all plastomes lacked ycf1 and ycf2 genes. Gaultheria plastomes shared substantially contracted SSC regions that contained only the ndhF gene. Moreover, plastomes of Gaultheria leucocarpa var. yunnanensis contained an inversion in the LSC region and an IR expansion to cover the ndhF gene. Multiple rearrangement events apparently occurred between the Gaultheria plastomes and those from several previously reported families in Ericales. Our phylogenetic reconstruction using 42 plastomes revealed well-supported relationships within all nine Gaultheria species. Additionally, seven mutational hotspot regions were identified as potential DNA barcodes for Chinese medicinal wintergreens. Our study is the first to generate complete plastomes and describe the structural variations of the complicated genus Gaultheria. In addition, our findings provide important resources for identification of Chinese medicinal wintergreens.

Chinese medicinal wintergreen plants differ in the type and content of chemical components, for instance, the content of wintergreen oil (Mukhopadhyay et al., 2016;Luo et al., 2021), but are not easily distinguished by morphological characteristics, particularly when processed into a pharmaceutical preparation. For example, one class of components that varies among Chinese medicinal wintergreens includes lignan glucosides, which have the highest content in G. leucocarpa var. yunnanensis and are known to possess strong anti-arthritic activity (Ma et al., 2001(Ma et al., , 2002Cheng et al., 2009;Hu et al., 2020). In addition, G. longibracteolata is used to treat pulmonary heart disease due to its unique phosphodiesterase-4 (PDE4). However, misidentification of G. longibracteolata for G. leucocarpa var. yunnanensis or G. fragrantissima has resulted in the misuse of G. longibracteolata to treat rheumatism and inflammation (Luo et al., 2021). This misuse of wintergreen plants may decrease medicinal effects or even cause medical malpractice. Therefore, accurate species identification of Chinese medicinal wintergreens is necessary.
The rapid development of high-throughput sequencing technology (Moore et al., 2006) has increased the availability of plastomes and lead to improved phylogenetic resolution at the interspecific level for many genera, e.g., Acacia, Atractylodes, Rhododendron and Salvia (Williams et al., 2016;Wang et al., 2021;Wu et al., 2021a;Fu et al., 2022). Zhang et al. (2017) resolved the inter-specific relationships of 19 species of the series Trichophyllae within Gaultheria, in which one medicinal wintergreen, G. sinensis, is phylogenetically placed. However, this study failed to obtain the complete plastid genome from long PCR sequences due to the presence of long repeats of a certain size, which hinder plastid genome assemblies of most species within Ericaceae (Fajardo et al., 2013;Fu et al., 2022). Structural variations in plastome, including IR expansion or contraction, inversion, and gene duplication or loss, together with DNA sequences can be used as effective tools for plant phylogenetic analyses (Ahmed et al., 2013;Daniell et al., 2016;Liu et al., 2020). Moreover, mutations in plastomes are usually concentrated in certain regions that form mutation 'hotspots' (Dong et al., 2012), suggesting that plastome sequences can be used to exploit short DNA barcodes and markers. To date, complete plastid genomes and plastome structural variation in Gaultheria have yet to be elucidated.
In this study, we sequenced the plastid genomes of 42 individuals from all nine recognized medicinal wintergreen species within Gaultheria. The aims of this study were to (1) resolve phylogenetic relationships among the nine Chinese medicinal wintergreen species; (2) investigate the structural variation of plastid genomes in Gaultheria at both the population and species levels; and (3) propose potential DNA barcodes for future wintergreen species identification.
Genomic DNA was extracted using the modified CTAB method (Doyle and Doyle, 1987), and the quality and quantity of genomic DNA were checked. The Illumina HiSeq 2500 platform was used to sequence the 500 bp insert-size libraries and to generate 2e3 Gb pair-end reads of 150 bp in length for each sample.

Assembly, annotation, and validation of complex repeat regions
To assemble the clean reads (processed reads) into complete plastomes, de novo assembly was performed with the GetOrganelle toolkit (Jin et al., 2020b). Connection and annotation were subsequently conducted using Bandage 0.8.1 (Wick et al., 2015) and Geneious 9.1.4 (Biomatters Ltd., Auckland, New Zealand) with Vaccinium macrocarpon Aiton. (GenBank accession number: NC019616) used as the reference. Notably, some plastid regions could not be disentangled in Bandage 0.8.1 (Wick et al., 2015) due to complex or long-size repeats (Fig. 1). These regions were subsequently verified by PCR and Sanger sequencing, using newly designed primers (Table S1).

Detection of SSRs and repeat sequences
REPuter software (Kurtz et al., 2001) was used to identify forward (F), reverse (R), palindrome (P), and complementary (C) repeats in the resultant plastomes with the setting of a minimum repeat size of 30 bp and 90% or greater sequence identity (Hamming Distance ¼ 3). Tandem Repeats Finder v.4.04 (Benson, 1999) was used to detect tandem repeats in the plastomes, with parameters set to 30 bp for the minimum repeat length, 0% for maximum mismatches and excluding repeats up to 10 bp longer than contained repeats. Simple sequence repeats (SSRs) were predicted using MISA (Beier et al., 2017) with the parameters of mono-, di-, tri-, tetra-, penta-and hexa-nucleotides set as 10, 5, 4, 3, 3, and 3, respectively. The maximum length of the sequence between two SSRs to register as a compound SSR was set to 0 bp.

Structural variation of Gaultheria plastomes
To investigate the structural characteristics of the Gaultheria plastomes under a broader scale, we downloaded representative plastomes from six related families within Ericales (i.e., Ericaceae, Lecythidaceae, Pentaphylacaceae, Primulaceae, Styracaceae and Theaceae; detailed information is provided in Table S2), and compared the structure of these plastomes with our newly generated plastomes using the MAUVE plug-in in Geneious (Darling et al., 2004).

Phylogenetic analyses
We downloaded three plastomes of the genus Vaccinium as outgroups (Vaccinium bracteatum GenBank No.: LC521967, V. uliginosum, LC521968 and V. vitis-idaea, LC521969). MAFFT v.7.310 (Katoh and Standley, 2013) was used to align the plastome sequences. Gene-partitioned model estimated using PartitionFinder2 (Lanfear et al., 2017) were used in the maximum-likelihood (ML) and Bayesian inference (BI) analyses. ML analysis was conducted using IQ-tree 1.6.12 (Nguyen et al., 2015) with support values estimated by 10,000 ultrafast bootstrap replicates. BI analysis was conducted using MrBayes 3.2.6 (Ronquist et al., 2012). Two runs with four chains were performed in parallel for ten million generations, with trees being sampled every 1000 generations.

DNA barcodes screening
The nucleotide diversity of 42 plastomes from nine species was calculated based on sliding window analysis using DnaSP v.5.10 (Librado and Rozas, 2009), with the window length set to 600 bp and step size set to 100 bp. Regions with high Pi values (> 0.2) were picked out initially and checked within Geneious v.8.1.7 (Kearse et al., 2012) with PCR success rates considered. The selected regions at this stage were subsequently used in phylogenetic tree reconstruction using RAxML (Stamatakis, 2014). The discriminability of DNA barcodes was evaluated by reconstructing phylogenetic trees. Species in monophyly with a support rate greater than 80% was considered to be identified successfully (Barrett et al., 2016).

Characteristics of plastomes of the Chinese medicinal wintergreens
The plastomes of all sampled Gaultheria share a typical quadripartite structure (Fig. 2), with GC content ranging from 36.5% to 36.7%. The plastomes ranged from 174,337 bp to 182,430 bp, with the LSC region ranging from 106,745 bp to 107,848 bp, the SSC region from 76 bp to 3642 bp, and the IR regions from 31,954 bp to 37,449 bp (Table 1). The total number of annotated genes was 137 in all samples of G. fragrantissima, G. hookeri, G. longibracteolata, G. nummularioides, G. pyrolifolia, G. semi-infera and two samples of G. griffithiana (LL50 and LL51), whereas it was 135 in all samples of G. sinensis, one sample of G. griffithiana (X19), and 138 in all samples of G. leucocarpa var. yunnanensis. Although the total number of chloroplast genes differed, the number of unique genes in each plastome was identical. Each plastome encoded 110 unique genes, including 76 protein-coding, 30 tRNA and 4 rRNA genes. All the Gaultheria plastomes have substantially enlarged IR regions, and the SSC region contains only one ndhF gene; moreover, the IR region of all samples of G. leucocarpa var. yunnanensis also contain the ndhF genes. We annotated two extra copies of the rpl23 gene (not identical in sequences) in the IR regions of all Gaultheria plastomes except the samples of G. sinensis and one sample of G. griffithiana (X19). We also annotated two extra identical copies of the adjacent trnfM and rps14 genes in the IR regions, and one extra copy of the trnfM gene in the LSC region (not identical in sequence) in all Gaultheria plastomes.

IR boundary shifts and genome comparison
We aligned all the resultant plastomes to investigate the expansion and contraction of IR boundaries. The plastomes showed high stability at the IR/LSC boundaries but were more varied at IR/ SSC boundaries. The psbA gene spanned the IRb/LSC junction, with 236 bp on the LSC side adjacent to the matK gene and 826 bp on the IRb side adjacent to the trnH gene. For the IRa/LSC junction, the trnV gene was located in the LSC region 51e63 bp away from the junction, whereas the trnH gene was located in the IRa region 1175e1179 bp away from the junction. The boundaries of IR/SSC regions can be classified into two types. Type I: the ndhF gene was located in the SSC region, 15e261 bp away from the IRb/LSC boundaries and 261e1296 bp away from IRa/SSC boundaries. This boundary type was observed in all sampled Gaultheria species except G. leucocarpa var. yunnanensis (Fig. 2). Type II: the ndhF genes were located in the IR regions due to IR expansion, and located 256e313 bp away from the IRb/SSC junction and 256e313 bp away from the IRa/SSC junction; no genes were present in the SSC region (Fig. 2). This boundary type was only observed in all samples of G. leucocarpa var. yunnanensis. Moreover, six out of nine samples of G. leucocarpa var. yunnanensis (incl., LL59, LL66, X26, X27, X28 and X29) had a unique inversion from the psbJ gene to the trnE-UUC gene (about 14 kb in length) in the LSC region, which resulted in the adjacencies of psbJ with petA and trnE-UUC with psbB (Fig. 2).
In addition, we compared the plastome structure of Gaultheria with plastomes from Ericaceae and five other families of Ericales (Lecythidaceae, Pentaphylacaceae, Primulaceae, Styracaceae, and Theaceae; Table S2). We found that the gene order of Gaultheria plastomes was collinear with Vaccinium (Ericaceae) but not collinear with that of Rhododendron (Ericaceae) (Fig. 3). The gene order of plastomes from five families of Ericales, excluding Ericaceae, are collinear (Fig. 3); however, the gene orders among Gaultheria, Rhododendron, and those five families can be explained by more than ten inversions (Fig. 3).

Phylogenetic inference
The resulting ML and BI tree topologies were highly similar to one another. Our phylogenetic analysis using plastomes in which one of the IR copies was removed revealed that all samples of each wintergreen species clustered into a monophyletic group with moderate supports for Gaultheria hookeri (BS ¼ 69, PP ¼ 0.99) and G. semi-infera (BS ¼ 87, PP ¼ 1.00), and strong supports for seven other species (BS ¼ 100, PP ¼ 1.00; Fig. 5). G. leucocarpa var. yunnanensis was sister to the remaining eight wintergreen species. G. sinensis was sister, in turn, to a grade consisting of G. griffithiana, G. pyrolifolia, G. longibracteolata, G. nummularioides, G. fragrantissima, G. semi-infera, and G. hookeri.

Gene content and structural variations of Gaultheria plastomes
The gene content of plastomes of land plants is generally considered to be conserved (Kumar et al., 2016), and the loss or duplication of different genes has been found in a variety of plants (Mohanta et al., 2020). It has been inferred that genes lost in plastomes may have moved to the nuclear genome or the function of these genes has been replaced by nuclear-encoded genes (Bryant et al., 2011;Savage et al., 2013;Mohanta et al., 2020). We identified several protein-coding gene losses and duplications in the Gaultheria plastomes. Both the ycf1 and ycf2 genes were absent in all studied Gaultheria samples. These two genes usually have elevated substitution rates and may have undergone pseudogenization in some land plant lineages (Oliver et al., 2010;Wolf et al., 2011). ycf1 has also been found missing in grasses and some parasitic plants, and independent losses of the ycf2 gene have occurred in multiple angiosperm groups, including Poaceae, Podostemaceae, Geraniaceae, and non-photosynthetic Ericaceae (Huang et al., 2010;Jin et al., 2020a). However, the mechanism for the loss of these two genes in the plastome is still not well understood.
In contrast to gene losses, gene duplications, attributed to repeat regions or IR expansions, were also detected in Gaultheria plastomes (incl., ndhF, rpl23, rps14, and trnfM). The ndhF gene was duplicated in all the samples of G. leucocarpa var. yunnanensis due to IR expansion. The chloroplast NAD(P)H dehydrogenase (ndh) complex is involved in photosystem I (PSI) cyclic electron transport and chlororespiration (Peng et al., 2011), and their mutation rates are high and sensitive to environmental conditions and stress (Silva et al., 2016). G. leucocarpa var. yunnanensis has a widespread Fig. 3. Syntenies and rearrangements detected in the complete plastomes of nine Chinese medicinal wintergreens (Gaultheria) and the genera used for comparison from six related families within Ericales. Boxes of the same color connected with lines indicate the local collinear blocks and represent syntenic regions. Sequence identity similarity profiles are represented by histograms within each block, and sequences transcribed in reverse directions are indicated by the colored blocks above and below the center line. The taxa studied are divided into six types of plastome arrangements: Type 1, samples LL59, LL66, X26, X27, X28, and X29 of G. leucocarpa var. yunnanensis; Type 2, samples LL61, LL63, and LL65 of G. leucocarpa var. yunnanensis; Type 3, all samples of G. fragrantissima, G. griffithiana, G. hookeri, G. longibracteolata, G. nummularioides, G. pyrolifolia, G. semi-infera, and G. sinensis; Type 4, Vaccinium from Ericaceae; Type 5, Rhododendron from Ericaceae; and Type 6, genera studied from Pentaphylacaceae, Styracaceae, Primulaceae, Lecythidaceae, and Theaceae.
Y.-L. Xu, H.-H. Shen, X.-Y. Du et al. Plant Diversity 44 (2022) 519e529 distribution in southern China, while other wintergreens inhabit the Himalaya-Hengduan regions. This difference in distribution implies that the duplication of ndhF gene is related to environmental adaptation in this species. In addition, plastomes with duplicated ndh genes usually exhibit extensive rearrangements and a higher frequency of pseudogenes (Martin et al., 1996;Casano et al., 2001;Paredes and Quiles, 2013). Consistent with this phenomenon, we identified an inversion in the LSC of some samples of G. leucocarpa var. yunnanensis. It is rare that duplication of plastid genes does not involve movement of IR boundaries (Raubeson and Jansen, 2005;Bock, 2007;Guisinger et al., 2011). In fact, only a few instances have been documented, for example, the duplication of clpP gene in Silene L. and Lychnis L. (Erixon and Oxelman, 2008), the duplication of psbJ and trnI-CAU in Trachelium Tourn. ex L. (Haberle et al., 2008), and the duplication of rpl23 and three tRNA genes in Jasminum L. (Lee et al., 2006). In this study, we identified duplication of three genes (i.e., rpl23, rps14 and trnfM) in Gaultheria plastomes that were independent of IR boundary shifts. The long repeats located in the LSC (between psaB and trnG) and both IR regions (between rrn16 and trnH) that covered rps14 and trnfM genes seemingly generated the three copies of these two genes. Additionally, we identified a fourth copy of trnfM in the LSC region adjacent to atpI, but with a slightly divergent sequence. By checking against other plastomes within Ericales, we inferred that the rps14 and trnfM genes between psaB and trnG are likely the original copies. Similarly, the rpl23 located in the LSC region (between trnI and rpl2) was not identical to the other two copies in the IR regions (between trnH and rps16); in this case, the LSC copy is likely the original copy.
Variation in plastome structure most typically involves expansion or contraction of the IRs into or out of adjacent single-copy regions (Ravi et al., 2008;Downie et al., 2015;Ye et al., 2018). The trnN-GUU gene is considered the ancestral IR/SSC endpoint that has been retained in most land plants, and the trnV-GAC gene represents the ancestral IR/LSC endpoint among most land plants (Zhu et al., 2016). Plastome IR expansions have been found in different plant groups, including Geraniaceae (Guisinger et al., 2011), Euphorbiaceae (Li et al., 2017), Solanaceae (Amiryousefi et al., 2018), Bignoniaceae (Thode and Lohmann, 2019) and Leguminosae (Bai et al., 2021). All plastomes of Gaultheria fragrantissima, G. griffithiana, G. hookeri, G. longibracteolata, G. nummularioides, G. pyrolifolia, G. semi-infera and G. sinensis expanded toward the SSC region from trnN to rpl32, and all samples of G. leucocarpa var. yunnanensis expanded toward the SSC region from trnN to ndhF. Different mechanisms have been proposed to explain IR expansions, such as gene conversion or double-strand DNA breaks (Goulding et al., 1996;Wang et al., 2008). In addition, some studies point out that the contraction and expansion of IR regions is a major contributor of plastomes size variation . For example, the IR regions of species from the genera Acorus and Magnolia are not expanded, and their plastome sizes are between 150 and 160 kb (Zhu et al., 2016), whereas the plastome size of Chinese medicinal wintergreens is between 176 and 182 kb. Our study validated that expansion of the IR region contributes to the increased genome size of Gaultheria species.
Similar to IR boundary shifts, inversions also represent essential mechanisms for plastome rearrangements, which contribute to the structural diversification of plant plastomes (Wicke et al., 2011;Jansen and Ruhlman, 2012). In nearly all photosynthetic land plants, plastomes are highly conserved in structural organization and gene arrangement (Jansen and Ruhlman, 2012;Zhu et al., 2016;Ye et al., 2018). Despite the commonly held view that plastomes have a conserved structure and sequence, they have been found to exhibit considerable variation in an increasing number of taxa Ye et al., 2018). In our study, we found that there was one inversion in the LSC regions of six samples (LL59, LL66, X26, X27, X28 and X29) within Gaultheria leucocarpa var. yunnanensis, which disrupts the canonical order of the plastome genes. In addition, except for the inversion in the above samples in G. leucocarpa var. yunnanensis, we found that there were apparent multiple rearrangement events between the Gaultheria plastomes and those from the other five families (Lecythidaceae, Pentaphylacaceae, Primulaceae, Styracaceae and Theaceae) in Ericales. Our study provides new evidence of the variation between Ericales plastomes.

Abundant and variable repeat sequences
The long repeats of Gaultheria plastomes were characterized in terms of the number, length and distribution. Previous studies indicated that long repeats mainly occur in non-coding DNA sequences (Raubeson and Jansen, 2005). However, we detected long repeats in coding regions of Gaultheria plastomes (e.g., infA, ndhA, ndhF, ndhI, petD, rpl22, rpl23, rpl32, rpoA, rpoC1, rps14, rps3, trnfM, and trnI). The most abundant sizes of long repeats among angiosperms are on average smaller than 50 bp (Raubeson and Jansen, 2005). With the increase in the number of plastomes sequenced, some longer repetitive sequences (approximately 100 bp) have also   6. Sliding window analysis of the complete plastomes of nine Chinese medicinal wintergreens (Gaultheria). X-axis: position of the window midpoint, Y-axis: nucleotide diversity within each window (window length: 600 bp, step size: 100 bp). Genic regions with Pi values > 0.2 are marked; region names with asterisk superscripts denote recommended DNA barcodes, while other regions are considered unsuitable for DNA barcodes due to the lack of either species discrimination or PCR amplification feasibility.
Y.-L. Xu, H.-H. Shen, X.-Y. Du et al. Plant Diversity 44 (2022) 519e529 been identified within some taxa. For example, the lengths of repeats range from 30 to 90 bp in Orchidaceae (Dong et al., 2018), 91 to 132 bp in Poaceae (Zhang et al., 2011), and 30 to 307 bp in Polypodiaceae . Generally, the number of repeats in the plastome is less than 50; for example, there are approximately 30 long repeats in Cymbidium (Yang et al., 2013), and 30 long repeats in Debregeasia (Wang et al., 2020). However, the longest long-repeat of each Gaultheria sample ranged from 706 bp to 1722 bp, and the number of long repeats was up to 259, which is significantly more abundant and variable than for many previously published plastomes.
Previous studies have shown that long repeats play a role in sequence rearrangement and variation in plastomes (Cavalier-Smith, 2002;Asano et al., 2004;Timme et al., 2007;Yang et al., 2013), and that large numbers of repeats can damage the stability of the plastome structure, for example, by producing inversions (Alexandre and Normand, 2010). Approximately 25 long repeats were located in the psbB-psbJ and trnE-petA regions of Gaultheria plastomes; thus, we speculate that this is one of the reasons for the 14-kb inversion in some samples of G. leucocarpa var. yunnanensis. Previous studies suggest that overall repeat content is positively correlated with the degree of genomic rearrangements (Wu et al., 2021b). This could be one of the reasons for the rearrangement of the plastomes in Ericaceae compared to other families within Ericales.

Phylogenetic relationships and delimitation of Chinese wintergreen species
Previous phylogenetic studies using a multiple-gene approach and international DNA barcoding (ITS, Waxy, Lfy, matK, rbcL, rpl16, trnH-psbA, trnL-trnF and trnS-trnG) failed to resolve the phylogenetic relationships among Chinese wintergreen species (Lu et al., 2010(Lu et al., , 2019bRen et al., 2011). In this study, using plastomes in which one of the IR copies had been removed to reduce redundancy, we resolved the phylogenetic relationships among all nine Chinese medicinal wintergreens.
Previous phylogenetic analysis recovered three major clades in Chinese Gaultheria, i.e., the Leucothoides clade, the Trichophyllae clade, and the Gymnobotrys clade, with the Leucothoides clade sister to the Trichophyllae clade (Lu et al., 2019b). In this and other analyses, G. leucocarpa var. yunnanensis, which diverged earlier than other studied species, was resolved as a member within the Gymnobotrys clade; G. fragrantissima, G. griffithiana, G. hookeri, G. longibracteolata, G. nummularioides, G. pyrolifolia and G. semiinfera were placed within the Leucothoides clade, but with poorly supported infraspecific relationships, and G. sinensis was placed in the Trichophyllae clade (Lu et al., 2010(Lu et al., , 2019aZhang et al., 2017). Our research supports this phylogeny, with G. sinensis sister to the grade comprising G. fragrantissima, G. semi-infera, G. griffithiana, G. hookeri, G. nummularioides and G. longibracteolata. We found that G. nummularioides had a sister relationship with the grade of G. hookeri, G. semi-infera and G. fragrantissima, and the clade comprising G. hookeri, G. semi-infera and G. fragrantissima was the most derived, which is consistent with the results of Lu et al. (2010). However, the relationship among the latter three species was poorly resolved in previous research (Lu et al., 2010;Fritsch et al., 2011), whereas our study suggests that each of the three species form a monophyletic lineage. In contrast to previous studies, we succeeded in clarifying the phylogenetic relationships among these nine species within the genus Gaultheria.

Candidate DNA barcodes
Chinese medicinal wintergreens are an important commodity in China, and the lack of genomic resources for Gaultheria has been the main obstacle to taxonomical and genetic analysis, as well as identification and conservation. Common barcodes are not effective for these species (Lu et al., 2010;Ren et al., 2011); although plastomes have been successfully used as a plant super-barcode to distinguish these species, plastome super-barcodes are costprohibitive and require complex bioinformatics processing (Shen et al., 2019). Short DNA barcodes, or markers, may provide an effective and economical alternative for identification of target plants (Chase and Fay, 2009;Chen et al., 2010;Shen et al., 2019). Fortunately, the mutation events in plastomes are not randomly distributed within the sequence and are concentrated in certain 'hotspot' regions; therefore, specific barcodes or markers can be exploited from the plastome (Dong et al., 2012). In this study, we proposed eleven regions with high nucleotide variability (Pi) values that might have a high potential for species delimitation in Chinese medicinal wintergreens. These mutational hotspots have the potential to resolve taxonomic issues in the genus and to be used in the future as barcodes for species identification. An ideal DNA barcode should be easily retrievable with a universal primer pair, an appropriately short sequence length to facilitate DNA extraction and amplification (Shen et al., 2019). The intergenic regions psbB-psbJ and trnfM-rrn16 possess high nucleotide variability (Pi) values, however, their PCR success rates were rather low (8.3% and 0 for each, respectively; tested with six samples of each wintergreen by using four primer pairs); thus, these two regions were excluded for the consideration of DNA barcodes. The rpoA-rpoB region was too long (approximately 6500 bp) for Sanger sequencing. Therefore, we finally selected seven regions (infA, ndhF, psaI-trnM, rps3-rpl16, trnQ-psaA, trnV-rps12 and ycf3-trnK) that have the potential to develop into DNA barcodes. Our results show that the trnV-rps12 region provided the best resolution for identifying all nine species of Chinese medicinal wintergreens. It was followed by trnQ-psaA region, which identified eight out of nine species. The infA gene provides the least resolution, which allowed the identification of only four out of nine species. In addition, eight divergence hotspot regions, including rpl36-infA, trnF (GAA)-ndhJ, trnD (GUC)-psbM, trnL (UAA)-trnF (GAA), trnT (UGU)-trnL (UAA), rpl23b, rps3b and petN-trnC (GCA), were recommended as candidate DNA barcodes for species identification of the Trichophyllae clade . In summary, different sets of DNA barcodes were detected for delimitation of species from the Trichophyllae clade and the Leucothoides clade.

Conclusions
Our research showed that the complete plastome could be used as a plant super-barcode to distinguish closely related medicinal wintergreens successfully, and that even specific molecular markers can be screened based on plastome sequences to distinguish the species. In addition, several unique characteristics of plastomes in Gaultheria, such as the loss of genes, rare gene duplication events, inversions (within G. leucocarpa var. yunnanensis), multiple rearrangement events, and remarkably long repeats, reveal new perspectives for understanding the dynamics of angiosperm plastomes. Furthermore, reticulate hybridization or species radiation that may have occurred during the evolutionary history of Gaultheria was not revealed by plastome sequence data, and further studies with increased taxa sampling and nuclear DNA sequences are necessary.

Author contributions
LL and XYD conceived the work. YLX performed the experiments and analyses. HHS contributed materials/analysis tools. YLX and XYD wrote the paper. XYD and LL revised the paper. All authors approved the final paper.

Declaration of interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.