Two Retrotransposon Elements in Intron of Porcine BMPR1B Is Associated with Phenotypic Variation

It has been established that through binding to bone morphogenetic proteins (BMPs), bone morphogenetic protein receptor I B (BMPR1B) can mediate transforming growth factor β (TGF-β) signal transduction, and is involved in the regulation of several biological processes, such as bone and muscle formation and homeostasis, as well as folliculogenesis. Also known as FecB, BMPR1B has been reported as the major gene for sheep prolificacy. A number of previous studies have analyzed the relationship between single nucleotide polymorphisms (SNPs) in this gene and its related performance. In recent years, with the illustration of the effect of retrotransposon insertion on the expression of the proximal genes or phenotypic variation, retrotransposon insertion polymorphisms (RIPs) have been used as a novel type of molecular marker in the evaluation of evolution, population structure and breeding of plant and domestic animals. In this study, the RIPs in porcine BMPR1B gene were excavated, and thereafter verified using a comparative genome and polymerase chain reaction (PCR). The potential effects of phenotype, gene expression and functions related to RIPs were also explored. The results showed that 13 distinct RIPs were identified in introns of porcine BMPR1B. Among these, only BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13 displayed a close relationship with the growth traits of Large White pigs. Moreover, the total number of BMPR1B-SINE+/+-RIP9 individuals born was found to be significantly higher than that of SINE−/− (p < 0.05). These two RIPs showed an obvious distribution pattern among Chinese indigenous breeds and Western commercial breeds. The expression of BMPR1B in ovaries of adult BMPR1B-SINE+/+-RIP9 Sushan pigs was found to be significantly higher in comparison to those of BMPR1B-SINE−/−-RIP9 (p < 0.05). SINE insertion of BMPR1B-SINE-RIP9 and LINE insertion of BMPR1B-LINE-RIP13 were observed to significantly increase the activity of Octamer binding transcription factor 4 (OCT4) minipromoter in CHO and C2C12 cells (p < 0.01). Therefore, these two RIPs could serve as useful molecular markers for modulating the growth or reproductive traits in assisted selection of pig breeding, while the mechanisms of the insertion function should be studied further.


Introduction
Bone morphogenetic protein receptor I B (BMPR1B) is a member of the bone morphogenic protein (BMP) receptor family, which can effectively bind BMPs to mediate transforming growth factor β (TGF-β)-induced signal transduction events [1], and is involved in the regulation of important biological processes, such as bone and muscle formation and homeostasis [2,3], folliculogenesis [4], etc. The protein encoded by BMPR1B gene contains four conserved residues, which can play different roles in defining BMP2, 6, 7 and GDF5 ligand affinity [5]. The variation and expression change of the BMPR1B have been related to the human ovarian insufficiency and defects of folliculogenesis [6,7], Pierre Robin sequence [8], isolated acromesomelic dysplasia [9], neuroblastoma [10], growth plate defects [11] and cancers [2].
For domestic animals, the point mutation of A746G in the coding region of sheep BMPR1B, also known as FecB, was reported to cause a significant increase in the ovulation rate and Booroola phenotype in 2001 [12][13][14]. Since then, a number of studies have been performed to decipher the relationship between this causal mutation and the prolificacy of sheep or goats [15][16][17]. The expression difference among these mutation allele individuals has been identified, in order to further elucidate the mechanism of the A746G mutation on biological processes [18][19][20]. The gene editing technology on this point mutation was used to successfully obtain the FecB-mutated sheep [21]. However, the transgenic pigs generated by introducing 746G in BMPR1B showed a normal litter size performance [22]. In addition, other point mutation and indels in the BMPR1B gene were found to significantly affect the litter size of sheep [23,24]. Overall, this gene has been thought as the major gene involved in regulating sheep prolificacy, and conformed connecting with pig prolificacy [25]. Up to now, there have only been few studies reported on the retrotransposon insertion polymorphisms (RIPs) of BMPR1B in pig breeding.
Retrotransposons, which are the dominant type of transposable elements (TEs), can effectively contribute to one-third to one-half of the mammalian genome [26]. They are an interesting genome component that can mobilize and generate new copies of themselves, as they are reversely transcribed into the genome using the pattern of "copy and paste". The cis-regulatory sequences contained in TEs could recruit the host-encoded factor, such as transcription factors (TFs), and ensure their amplification within the new loci of host [27]. Therefore, TE insertion can affect the proximal genes by functioning as enhancers, promoters, silencers and boundary elements contributing to the genome architecture [28,29].
According to the presence of long terminal repeats (LTRs) or non-LTR, retrotransposons can be further classified into two types. Endogenous retroviruses (ERVs) are the main part of LTRs. Non-LTR retrotransposons have been divided into long interspersed elements (LINEs) and short-interspersed elements (SINEs). Based on the retrotransposon insertion polymorphism, several types of molecular markers have been developed in plant breeding [30,31]. Recently, RIPs were also identified as a causative mutation of the phenotype, or closely related to the performance of domestic animals [32][33][34][35][36][37][38][39]. Considering the importance of the BMPR1B protein in several biological process, the RIPs in BMPR1B were excavated and verified. The potential relationship of these RIPs with pig growth, as well as the reproductive performance and the mechanism of these insertions, were identified to find the effective and novel molecular markers that could be used in pig breeding.

Animals and DNA Isolation
Ear or blood samples were collected from the twelve pig breeds or population including Duroc, Landrace, Large White, Sujiang, Diannan small ear, Erhualian, Wuzhishan, Bama, Tibetan, Meishan, Fengjing and wild boars. Genomic DNA was extracted using a standard phenol-chloroform extraction procedure and stored at −20 • C. DNA pool from 12 breeds or population were prepared for RIPs verification, with 6 individuals for each breed. A certain number of breeds, including Duroc, Landrace, Large White, Erhualian, Wuzhishan, Sushan and Sujiang were prepared for analyzing the population diversity. The correlation analysis was performed between the different RIP genotypes and the performance records of thickness of backfat, age at 30 kg bodyweight, age at 100 kg bodyweight, first litter size, number born alive, litter weight of Large White pigs. Both the origin of breeds and number have been shown in Table S1. Among them, Large White, Duroc and Landrace are three lean type commercial breeds. Fengjing, Meishan, Erhualian and Wuzhishan are four fat type indigenous breeds. Sujiang is a cross breed which originated from Duroc, Jiangquhai and Fengjing. Sushan is also a cross breed developed from Meishan, Erhualian, Landrace and Large White in Jiangsu Province [40].

RIPs Excavation
The sequence of BMPR1B (Gene ID: 396691) with the flanking region (5 kb 5 upstream and 3 kb 3 downstream) were downloaded from NCBI (https://ncbi.nlm.nih.gov/, accessed on 17 May 2021). As the bait, BMPR1B sequences from the fourteen assembled non-reference genomes (Bama, Pietrain, Jinhua, cross-breed of Yorkshire/Landrace/Duroc, Landrace, Hampshire, Wuzhishan, Berkshire, Large White, Tibetan, Meishan, Bamei, Gottingen and Rongchang) were downloaded from the NCBI database (https://blast.ncbi.nlm. nih.gov/Blast.cgi, accessed on 19 May 2021). All the sequences were aligned to excavate structural variations using the ClustalX (Version 2.0). In addition, large structural variations (more than 50 bp and presenting more than one population) were retained for the retrotransposon annotation using RepeatMasker (Versions 4.0.7). The structural variations with over 60% similarity to the sequences in RepeatMasker were considered as the predicted retrotransposons and used for further identification by PCR.

RIPs Identification by PCR
The primers spanning the predicted retrotransposons were designed using Oligo7. The primer sequences and annealing temperature are shown in Table S2. All the primers in this study were synthesized by Beijing Liuhe Huada Gene Technology Co. Ltd. (Beijing, China). PCR was carried out using ABI PRISM 3100 Genetic Analyzer (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA), and the reaction system was prepared as described previously by using the 12 DNA pools as template [38].

Population Diversity Based on RIPs
BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13 were used to illustrate the genetic diversity of the represented breeds, including the lean-type commercial breeds (Large white, Landrace and Duroc), fat-type indigenous breeds (Fengjing, Meishan, Erhualian and Wuzhishan) and Jiangsu cross breeds (Sujiang and Sushan). Genotype frequency, allele frequency and Hardy-Weinberg test of RIPs were analyzed by POPGENE 1.32 [41]. Polymorphic information content (PIC) was calculated by the following formula: The primers (Table S2) for the pig BMPR1B (mRNA XM_021100423.1) and GAPDH (mRNA XM_003126531.5) were designed based on the sequences available in GenBank Using Oligo7 (accessed on 3 June 2021).
Four Sushan individuals for each genotype of BMPR1B-SINE-RIP9 and for LINE +/+ and LINE +/− individuals of BMPR1B-LINE-RIP13 were slaughtered to obtain the ovary, backfat and longissimus dorsi. Thereafter, the tissue RNA was extracted using TRIzol (Tiangen, Beijing, China), and the quantity and quality of total RNA were monitored using a spectrophotometer (NanoPhotometer N60 Touch, Implen Gmbh, Munich, Germany). The cDNA was prepared using Takara Kit (Takara, Dalian, China). A real-time quantitative polymerase chain reaction (qPCR) was performed in a 25µL reaction mixture, containing 12.5 µL TB Green ® premix Ex TaqTM II (2×, Biomedical, Dalian, China), 1 µL each forward and reverse primer (10 ng/µL), 8.5 µL ddH2O, and 2 µL cDNA (100 ng/µL). The reaction conditions used were as follows: pre-denaturation at 95 • C for 30 s; denaturation at 95 • C for 3 s; and extension at 60 • C for 30 s (40 cycles). The relative expression levels of BMPR1B mRNA were analyzed using the 2 −∆∆Ct method, and the GAPDH gene was used as the internal reference.

Dual-Luciferase Reporter Assay
The insertion fragments with the double restriction enzyme cutting site of BMPR1B-SINE-RIP9 (273 bp) and BMPR1B-LINE-RIP13 (804 bp) were cloned from the Sushan genomic DNA (the primers have been listed in Table S2), and verified by sequencing. The fragments were then inserted into the pGL3-basic vectors (Promega, Madison, WI, USA). Furthermore, Oct4 mini promoters were cloned from the pTol2-Oct4-mCherry vector [42] and inserted into pGL3-basic vectors with or without the insertion fragment to evaluate the potential enhancer or depressor activity. The correct recombinant plasmids were designated as pGL3-SINE-Oct4-basic and pGL3-LINE-Oct4-basic. Using a Lipofectamine 3000 reagent (Invitrogen, Carlsbad, CA, USA), 2 × 10 4 CHO and C2C12 cells were plated and transfected with the plasmids. The cells were collected after 48 h to evaluate the luciferase activity using the dual luciferase reporter system (Promega, Madison, WI, USA) according to the manufacturer's instructions. The luciferase activity was calculated by the following formula: Luciferase activity = Firefly luciferase activity Renilla luciferase activity

Statistical Analysis
All the data has been shown as mean ± standard error, and p < 0.05 was considered a significant difference. The results were processed by SPSS17.0 (SPSS, Chicago, IL, USA) using one-way analysis of variance with Duncan's test. The association analysis was also calculated by using one-way ANOVA in the SPSS software. The general linear model was used to determine the possible association between RIPs and growth and reproductive traits. Equation 1 for growth trait and 2 for reproductive traits used were as follows: (1) y ijk = u + g i + a j + s k + e ijk , (2) y ijk = u + g i + s k + e ik , where y ijk was the phenotypic observation, u was the mean, g i was the fixed effect of genotype, a j was the fixed effect of gender, s k was the fixed effect of born weight, and e ijk or e ik was the residual effect.

Animal Welfare
All the treatments and protocols involving animals in this study were strictly performed in accordance with the guidelines of the Animal Experiment Ethics Committee of Yangzhou University (approval number: YZUDWSY2018-12).

Thirteen RIPs Were Verified by PCR in the Pig BMPR1B Gene
By RepeatMasker annotation, 38 SINE, 17 LINE and 6 LTR insertion polymorphisms were predicted in the BMPR1B gene among the reference genome, and another 15 nonreference genomes. The PCR was employed to identify these predicted RIPs by using the specific primers (Table S1), and 13 RIPs were verified in DNA pool ( Figure 1A). They were all distributed in the introns of BMRP1B gene ( Figure 1B). Most of the insertion loci were relatively not conserved ( Figure 1C). The details, including the insertion sequence size, direction and subfamily of retrotransposon, have been listed in Table 1. These RIPs were named as BMPR1B- and BMPR1B-LINE-RIP13, which corresponded to labels 1-13 in Figure 1B.  Note: "+" indicates forward insertion, "−" indicates reverse insertion.

Association between RIPs and Phenotype of Large White Pigs
The association of 13 RIPs with the phenotype of Large White pigs was analyzed, and BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13 were found to be closely associated with the growth and reproductive traits of Large White pigs. As shown in Table 2, the age at 100 kg body weight of SINE −/− individuals for BMPR1B-SINE-RIP9 was observed to be significantly lower than those of SINE +/− individuals (p < 0.01). The age at 30 kg body weight and 100 kg body weight of LINE -/individuals for BMPR1B-LINE-RIP13 was noted to be significantly lower than those of LINE +/+ individuals (p < 0.05). The first litter size of SINE +/+ individuals for BMPR1B-SINE-RIP9 of Large White pigs was significantly higher as compared to those of SINE −/− individuals (p < 0.05) ( Table 3). Note: +/+: homozygote with retrotransposon insertion; +/−: heterozygote with retrotransposon insertion; −/−: homozygote without retrotransposon insertion. The same superscript letters in the same column (a and ab or b and ab) mean that the difference between groups was not significant. Different superscript lowercase letters (a and b) indicated a significant difference between groups (p < 0.05). Different superscript capital letters (A and b) indicated an extremely significant difference between groups (p < 0.01). Note: +/+: homozygote with retrotransposon insertion; +/−: heterozygote with retrotransposon insertion; −/−: homozygote without retrotransposon insertion. The same superscript letters in the same column (a and ab or b and ab) mean that the difference between groups was not significant. Different superscript lowercase letters (a and b) indicated a significant difference between groups (p < 0.05).

BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13 Distribution in the Different Pig Breeds
The potential distribution and genetic parameters of BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13 in lean-type commercial breeds (Large White, Landrace and Duroc), fat-type indigenous breeds (Fengjing, Erhualian, Meishan, Jinhua and Wuzhishan) and cross breeds in Jiangsu province (Sujiang and Sushan) are shown in Table 4. The SINE + frequency of BMPR1B-SINE-RIP9 was dominant in the lean-type breeds (Large White and Landrace) and cross breeds (Sushan and Sujiang) while the SINE − was dominant in fat-type breeds (Wuzhishan). The LINE + frequency of BMPR1B-LINE-RIP13 was obviously found to be higher in lean-type breeds (Duroc, Large White and Landrace) in comparison to cross breeds (Sushan and Sujiang) and fat-type breeds (Jinhua and Erhualian). Except for Large White pigs, all RIPs in the different breeds were in Hardy-Weinberg equilibrium. BMPR1B-SINE-RIP9 in fat-type indigenous breeds (Wuzhishan) and Sujiang were low polymorphic, while lean-type breeds (Duroc, Large White and Landrace) and Sushan were medium polymorphic. BMPR1B-LINE-RIP13 in the cross pigs (Sushan and Sujiang) and fat-type breeds (Jinhua and Erhualian) were low polymorphic, but in lean-type breeds (Duroc, Large White and Landrace) they were medium polymorphic.

The Expression Pattern of Pig BMPR1B in Sushan Pigs for Two RIPs
The expression in ovary, longissimus dorsi and backfat of the different genotype individuals of six-month old Sushan pigs for BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13 was shown in Figure 2. For BMPR1B-SINE-RIP9, there was significant difference in the expression of BMPR1B among SINE +/+ and SINE −/− individuals in ovaries (p < 0.05). Two genotypes, LINE +/− and LINE −/− , were found in Sushan pigs for BMPR1B-LINE-RIP13. There was no significant difference noted between these two genotype individuals for the expression of BMPR1B in the three tissues of Sushan pigs (p > 0.05).  The values are shown as mean ± SD, and * showed p < 0.05. "+/+" means Homozygous insert. "+/−" means heterozygous insert. "−/−" means no insert. The values are shown as mean ± SD, and * showed p < 0.05.

SINE and LINE Insertion Might Act as an Enhancer in CHO and C2C12 Cells
The insertion fragments, which were 273 bp and 692 bp for BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13, respectively, were amplified from Large White genomic DNA and then inserted into the pGL3-Oct4-basic vector. The recombinant plasmids were named as BMPR1B-SINE9-Oct4-Luc + and BMPR1B-LINE13-Oct4-Luc + ( Figure 3C). These two plasmids were transfected into CHO and C2C12 cells, and the pGL3-Oct4-basic plasmid was used as a negative control ( Figure 3D). The luciferase activity of BMPR1B-SINE9-Oct4-Luc + and BMPR1B-LINE13-Oct4-Luc + was found to be significantly higher than that of pGL3-Oct4-basic negative control (p < 0.01) in CHO and C2C12, Therefore, the 273 bp fragment and 692 bp fragment might act as an enhancer to affect the Oct4 promoter activity in CHO and C2C12 cells.

Discussion
BMPR1B is a type 1B receptor and a pivotal member of the BMP receptor family. The protein encoded by the BMPR1B gene mainly participates in the regulation of the BMP signaling pathway. BMPs are secreted ligands in the TGFβ superfamily, and can perform various important roles, including neurogenesis, organogenesis, pregnancy, cancer, etc [43]. All porcine BMP genes were screened to find out useful RIPs, but nothing was found within the gene and upstream or downstream flanking in the previous reports. The similar results were also observed in GH, Leptin gene, which is necessary for maintaining the different physiological processes including growth, development and reproduction. The studies have shown that 85-90% of mouse and human protein coding genes contain transposon elements sequences in their introns [44]. Thus, it is possible that essential genes might have developed measures to control the retrotransposon mobility in order to ensure proper physiological function [45]. However, their receptors could present some RIPs, which lead to the sequence diversification of organisms. All identified RIPs were found in the intron of the porcine BMPR1B gene. In our previous study, only 0.63% of retrotransposons were inserted into the exon of protein coding genes, while 35.1% retrotransposons were inserted into the intron of protein coding genes [46]. The considerable transposon elements' under-representation in the coding regions could be attributed to a strong negative selection to avoid repercussions in the protein structure [47].
A number of studies have reported that BMPR1B is the major gene for determining the litter size of sheep. As a new-type molecular marker, 13 RIPs were used to analyze the correlation with the growth and reproductive traits of Large White pigs in this study. It was found that only BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13 were closely related to the growth traits of Large White pigs, while BMPR1B-SINE-RIP9 was closely linked to the first litter size of Large White pigs. Using high throughput sequencing technology, it has been reported that the selective sweep that included BMPR1B (chr8:134.05-134.15 Mb and chr8: 123.3-147.4 Mb), which overlapped with the QTL to BMPR1B, significantly affected the number of piglets born [25,48]. BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13 were within or very close to the selective sweep for litter size on the Chr8. Thereafter, the distribution of these two RIPs were evaluated in the representative pigs of three different type pigs, including fat-type, lean-type and cross breeds. Chinese indigenous pigs are famous for their high prolificacy, good meat quality, slow growth rate and low lean percentage in comparison with commercial pigs originating from Western nations [49][50][51][52][53]. There were obvious different distribution patterns noticed among these three types of pigs. Some nucleotide variants in the intragenic and intergenic regions indicated population characteristics among the different originated pig breeds [54][55][56][57]. BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13 showed the obvious variable distribution pattern among the different type of breeds, which may be attributed to the difference of artificial selection concerned with these pigs. Therefore, we focused our attention on these two RIPs.
BMPR1B defects in mice have been shown to cause infertility [58], and high expression of BMPR1B in ovaries has been related to highly prolific sheep and goats [20,59]. In this study, the SINE-insertion individuals exhibited higher expression of BMPR1B in ovary and greater first litter size in comparison to individuals without SINE insertion. The function of BMPR1B was diverse, and the protein contributed to the growth and development of the mammal. The two RIPs were also found to regulate the growth rate of Large White pigs. However, the expression difference was not found in back fat and muscle among these genotypes. It is possible that the expression of BMPR1B in the tissues of brain was relatively high compared with other tissues, including the ovary [20,23], which made it play a major role through modulating the hypothalamic-pituitary axis to affect the growth and development.
To explore the potential mechanism of retrotransposon insertion affecting the expression of BMPR1B, the dual-luciferase reporter assay was employed to evaluate the potential activities of insertion sequences. Retrotransposon can disperse vast amounts of promoters, enhancers, and even repressive elements [60,61]. In this study, the SINE and LINE insertion significantly increased the dual-luciferase activity in CHO and C2C12 cells. Since SINE insert is present upstream of CDS domain, it could act as an enhancer to increase the expression of BMPR1B in ovaries and affect the reproductive traits of Large White pigs. However, the mechanisms through which BMPR1B-LINE-RIP13 can affect growth traits is still unclear, but the function of LINE insertion should be further explored. Overall, combined analysis of the relationship between these two RIPs and related traits of Large White pigs could lead to useful molecular markers for assisted selection in pig breeding.

Conclusions
Thirteen distinct RIPs were identified in introns of porcine BMPR1B by comparative genome and PCR. Only BMPR1B-SINE-RIP9 and BMPR1B-LINE-RIP13 displayed a close relationship with the growth rate or first litter size of Large White pigs. These two RIPs exhibited an obvious distribution pattern that the SINE + frequency of BMPR1B-SINE-RIP9 and the LINE + frequency of BMPR1B-LINE-RIP13 was dominant in the commercial breeds in comparison to Chinese indigenous breeds, which may be attributed to the difference of artificial selection concerned with these pigs. The expression of BMPR1B in ovaries was altered significantly by the SINE insertion of BMPR1B-SINE-RIP9 in adult Sushan pigs (p < 0.05). SINE insertion of BMPR1B-SINE-RIP9 and LINE insertion of BMPR1B-LINE-RIP13 significantly increased the activity of OCT4 minipromoter in both CHO and C2C12 cells (p < 0.01). Therefore, these two RIPs could serve as useful molecular markers for predicting the growth or reproductive traits in the assisted selection of pig breeding.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/life12101650/s1, Table S1. Number and origin of pig breeds for BMPR1B RIPs detection. Table S2. The primers for PCR, vectors construction and q-PCR.