Deciphering the lncRNA and mRNA profiles of Min pig backfat after acute cold stress

ABSTRACT

Cold stress can promote fat metabolism and reduce fat deposition. The present study compared differentially expressed (DE) lncRNAs and mRNAs in the subcutaneous back fat of Min pigs subjected to 24 h acute cold stress. Transcriptome sequencing was used to identify DE lncRNAs and mRNAs. A lncRNA mRNA co expression network was constructed using the confirmed lncRNAs. GO and KEGG analyses were performed to investigate the biological processes involved in the response to acute cold stress. A total of 166 DE lncRNAs and 269 DE mRNAs were found. Several genes such as ANGPTL4, KLF11 and NUPR1 were related to lipid metabolism. Interferon stimulating genes (ISGS) such as RSAD2, ISG20, ISG15, DDX60, MX2, MX1 and IFI6 were all significantly up-regulated. In addition, genes related to blood flow, respiration, and oxidative stress were also up-regulated. Many terms related to defense response to a virus were found in the GO annotation. A lncRNA mRNA co expression network comprising 46 specific mRNAs and the top 10 up-regulated DE lncRNAs was constructed. According to the results, the interaction between lncRNAs and mRNAs was complex. This study provides genome-wide gene and lncRNA expression profiles in adipose tissue of pigs after acute cold stress.

KEYWORDS:

• Min pig
• adipose tissue
• transcriptome sequencing
• lncRNA
• acute cold stress

Introduction

With rapid economic development and the associated changes in lifestyle, obesity has become a global public health problem (Christensen 2020). Obesity results from an imbalance between caloric intake and energy expenditure, which promotes the expansion of adipose tissue, increases fat content, and increases the incidence of metabolic diseases such as type 2 diabetes. Obesity also contributes to cardiovascular and liver diseases (Marcelin et al. 2019). Browning of white adipose tissue (WAT) has been shown to reduce insulin resistance and the chronic inflammation caused by obesity (Shen et al. 2020; Kuryłowicz and Puzianowska-Kuźnicka 2020). Factors such as exercise and a cold environment can stimulate the browning of white adipose tissue (Piao et al. 2018; Martin et al. 2020). Activation of brown adipose tissue (BAT) and browning of WAT, although beneficial for the treatment of metabolic diseases, may also promote the disease progression in patients with malignant tumors (Dalal 2019).

LncRNAs are a form of non-coding RNA that play precise regulatory roles in development and gene expression. Different lncRNAs are involved in both the early and late stages of adipogenesis and are associated with complications of obesity. LncRNAs are gradually becoming known as regulatory factors (Rey et al. 2021). LncRNA Blnc1 forms a nucleoprotein transcription complex with Zbtb7b to promote the development and metabolic thermogenesis of brown and beige fat (Li et al. 2017). Lnc-BATE1 binds specifically to nuclear matrix factor hnRNPU and determines the differentiation of brown adipocytes (Hacisuleyman et al. 2014). LncRNA ORA promotes adipogenesis through the PI3K/AKT/mTOR signaling pathway (Cai et al. 2019). Alvarez-Dominguez et al. found 127 specific lncRNAs in brown adipocytes, many of which were related to key pathways that determine adipocyte differentiation, for example PPARγ and C/EBPα (Alvarez-Dominguez et al. 2015).

Pigs are both an economically important domestic animal and an emerging model in which many human genetic mechanisms and disease treatments can be analyzed (Pabst 2020; Zuidema and Sutovsky 2020). However, pigs differ significantly from humans and other mammals in their response to low temperatures. The Uncoupling protein 1 (UCP1) gene, a key gene in brown adipocytes production, has mutated over the course of evolution for reasons unknown. The mutation prevents the production of brown adipocytes (Berg et al. 2006). Therefore, how pigs maintain a constant body temperature at low temperatures is not completely understood. Previous studies have found that UCP3, a member of the same protein family as UCP1, can replace UCP1 in cold environments and promote white adipose tissue metabolism and heat production in pigs (Lin et al. 2017). If mouse UCP1 is overexpressed with Peroxisome Proliferator-activated Receptor γ coactivator-1 α (PGC-1α) in pigs, the white precursor adipocytes of pigs can be browned (Hou et al. 2018). However, to date, the gene expression changes in adipose tissue of pigs under cold environmental conditions have not been thoroughly examined.

Min pigs, an indigenous pig species in northern China, were used as experimental subjects. It is not only a good model to explore the genetics underlying economic traits (Wang et al. 2014), but also a good model of studying cold stress. The Min pigs belong to the fatty type and have innate cold tolerance, as they have lived in cold areas for a long time. Under the same cold environment, phenotypic and physiological changes in Min pigs have been more adaptive than those in Large White pigs (Liu et al. 2017), and thus the system provides an opportunity for a comparative study of adipose tissue metabolism under cold stress. In this study, transcriptional sequencing technology was used to analyze the expression profiles of genes and lncRNAs in the subcutaneous fat of Min pigs after 24 h cold stress.

Materials and methods

Experimental design and sample collection

Six six-month-old female Min pigs with an average weight of 60 ± 3 kg were provided by the Min pig breeding farm of Heilongjiang Academy of Agricultural Sciences. The six Min pigs had no genetic relationship and were randomly divided into two groups. One group was raised in a heat preservation house where the temperature was controlled at 20 ± 3°C (the control group). The other group was raised in a low temperature house where the temperature was controlled at −15 ± 3°C (the cold stress group). During the experiment, all pigs were given free access to food and water. After 24 h of cold treatment, all experimental pigs were slaughtered. The subcutaneous adipose tissue of the back of each pig was quickly removed and stored in liquid nitrogen.

Extraction and detection of total RNA from subcutaneous fat

Total RNA was extracted from the subcutaneous adipose tissue of Min pigs by the phenol/chloroform method. The purity of total RNA was detected by NanoDrop 2000 & 8000 micro spectrophotometers. The concentration and integrity of total RNA were detected by an Agilent 2100 Bioanalyzer and an Agilent RNA 6000 Nano kit.

Construction and sequencing of the lncRNA library

The lncRNA library was constructed with 3 μg total RNA from each sample. A Ribo-Zero™ Magnetic Gold kit was used to remove rRNA, and the fragmentation buffer was added to the reaction system to cut the RNA into short fragments. Using the fragmented RNA as a template, six-base random primers were used to synthesize the first strand of cDNA. The second strand of cDNA was synthesized by adding buffer solution, dNTPs, RNase H, and DNA Polymerase I. The cDNA was purified using a Qiaquick PCR kit. The cDNA fragments were end-repaired, after which base A and sequencing adaptors were added. The second strand of cDNA was digested by UNG enzyme, amplified by PCR, and purified. Finally, an Illumina HiSeq 2000 platform and paired-end sequencing strategy were used.

LncRNA analysis

The basic screening conditions for lncRNA in this study were as follows: (1) the length of a transcript being greater than or equal to 200 bp, and the number of exons being greater than or equal to 2. (2) The reads coverage of each transcript was calculated, and the transcripts with reads coverage less than 5 in all samples were screened out. (3) GffCompare was used to make comparisons with pig genome annotation files to screen out known mRNAs and other noncoding RNAs (i.e. rRNA, tRNA, snoRNA, or snRNA). (4) Potential lncRNAs were screened according to the class_code information from the comparison results (if class_code = ‘u’, this indicates lincRNA; if class_code =‘I’, this indicates intronic lncRNA; class_code =  ‘x’ indicates anti-sense lncRNA). After screening, four coding potential analysis software programs (CNCI, CPC, PFAM and CPAT) were used to judge whether the screened lncRNA had coding potential.

Analysis of differentially expressed genes

Gene expression levels were estimated using the FPKM method. The calculation formula was as follows: $\mathrm{F}\mathrm{P}\mathrm{K}\mathrm{M}\left(A\right)={\displaystyle \frac{{10}^3\times \mathrm{F}}{\mathrm{N}\mathrm{L}/{10}^6}}$, where F is the number of fragments uniquely aligned to gene A; N is the total number of fragments uniquely aligned to the reference gene, and L is the length of the exon region of gene A. The FPKM method can eliminate the influence of gene length and sequencing amount in calculating gene expression levels..

The genes with |log₂ Ratio| ≥1 and q < 0.05 were selected as having significant differential expression, and the numbers of DE genes in experimental groups were obtained. The DE genes in different groups were analyzed by a hierarchical cluster analysis using the R software (V3.1.1). The DE genes were annotated with NCBI, UniProt, GO, and KEGG databases to obtain the detailed description information. The lncRNA was annotated using ENSEMBL, NCBI, GENCODE and HGCN databases.

GO and pathway analyses

The online IGV software (http://www.broadinstitute.org/igv/download) was employed for DE genes’ GO annotation, with the categories being biological process (BP), cellular component (CC), and molecular function (MF). The redundant GO terms were removed by Revigo. The KEGG Orthology Based Annotation System (KOBAS) software (http://www.genome.jp/kegg) was employed for KEGG pathway enrichment analysis. P < 0.05 was considered as signifying that a function was significantly enriched.

Analysis of the lncRNA-mRNA co-expression network

The target genes regulated by lncRNAs through cis or trans modes were predicted. For those lncRNAs acting in cis mode, the genes at the adjacent positions (50 kb upstream and downstream) were screened as target genes. Those predicted to act in trans mode were screened according to the correlation coefficient between lncRNA and mRNA expression (correlation coefficient ≥ 0.9). The target genes of the top 10 up-regulated differential lncRNAs were analyzed. According to the identified target genes, a co-expression network diagram of lncRNA and mRNA was drawn using Cytoscape software (Shannon et al. 2003).

Quantitative reverse transcription- polymerase chain reaction (qRT-PCR)

Total RNA of every sample was isolated with TRIzol® reagent and reverse transcribed into cDNA using the PrimeScript RT reagent Kit gDNAEraser. In order to verify the accuracy of transcriptome sequencing data, 10 highly expressed lncRNAs were selected and verified by qRT-PCR. The primers used for qRT-PCR are listed in Table 1. qRT-PCR was performed using SYBR-Green quantitative PCR and the Roche LightCycler480 II system. The β-actin gene of pigs was used as an internal control. The PCR reaction system comprised 20 μL: SYBR, 10 μL; up and down primers (10 μmol/L), 0.6 μL; cDNA template, 1μL; RNase-free ddH₂O, 7.5 μL; Rox, 0.3 μL. The PCR reaction program was 95°C for 30 sec, 40 cycles at 95°C for 5 sec, and 60°C for 30 sec. The relative expression level of each lncRNA was calculated using the 2^−ΔΔCt method (Schmittgen and Livak 2008) and then statistically analyzed.

Table 1. Primers used for quantitative reverse transcription PCR.

lncRNA	Forward primer (5'−3’)	Reverse primer (5'−3’)
MSTRG.103678	TATCTTATGGACCCTGGTGAGC	CAATCAACTCTACCCAACTCCC
MSTRG.77107	ATTGAATCCGTTGATGGCAGTT	AGCGATGGAGGTGAACAATGAG
MSTRG.31867	CCCTCAAGCAACATAGATGGCAGAT	CTCCCGCCTCAGTAGCCGTTTT
MSTRG.85551	ACCACCATCCACTCATTTCCTT	CTCGCTCCAGCCTATGTTTCAG
MSTRG.107837	ACGGTTTAGTGGTTGGGTGAAG	GGAAGGGAACAATAGCCAGAAT
MSTRG.76112	AGCGGTAGAATGCTAATGGGTT	GGTTGCTCTTCCTGGAATGACT
MSTRG.82028	CTGGCAACGATGGGAGAAGGGAT	CAGGCTGAGTGGCTTGGCTAGGTC
MSTRG.44325	AGGGCGTTCCGAGAAGGTGCTG	CGTGTTTGTCCGTTACCGATGT
MSTRG.44154	GGGGATTACAATAGGGTAGAGC	TGGAAGAGGTTAGGGAACAGAG
MSTRG.146756	GCTGGGTCTGTAGTGGAGGAAT	ATGCTGGCAGGTCTGTAGGAAA
β-actin	CGGGACCTGACCGACTACCT	GGGCCGTGATCTCCTTCTG

Results

Analysis of lncRNAs and mRNAs in subcutaneous adipose tissue of Min pigs after acute cold stress

To obtain the expression profiles of lncRNAs and mRNAs in subcutaneous adipose tissue of Min pigs after acute cold stress, subcutaneous adipose tissue samples from three pigs under normal temperatures and three pigs subjected to acute cold stress were examined. The results revealed 7856 lncRNAs and 18556 mRNAs in total. Compared with the control group, 78 lncRNAs (yellow dots in Figure 1(a)) and 141 mRNAs (yellow dots in Figure 1(b)) were up-regulated (P < 0.05; log₂ Ratio > 1) in the cold stress group. In contrast, 88 lncRNAs (green dots in Figure 1(a)) and 128 mRNAs (green dots in Figure 1(b)) were significantly down-regulated (P < 0.05; log₂ Ratio < −1). All DE lncRNAs and mRNAs are presented in the hierarchical cluster analysis heat map (Figure 2(a,b), Supplement Table 1 and Table 2), in which ‘C group’ represents three subcutaneous adipose tissue samples from the acute cold stress group and ‘N’ represents three subcutaneous adipose tissue samples from the control group. Data can be obtained in NCBI database (accession No. PRJNA766153).

Figure 1. Volcano plot of DE lncRNAs and DE mRNAs. (A) The expression levels of mRNAs were compared between the acute cold stress group and the control group. (B) The expression levels of lncRNAs were compared between the acute cold stress group and the control group. The yellow dots represent transcripts significantly up-regulated in the acute cold stress group, and the green dots represent transcripts significantly down-regulated in the acute cold stress group based on the criteria (|log₂ Ratio| ≥ 1 and P < 0.05). The horizontal axis represents the fold change in expression (on a log2 scale), and the vertical axis represents the P-value (on a negative log10 scale).

Figure 2. Hierarchical clustering analysis of DE lncRNAs and DE mRNAs. (A) Hierarchical clustering analysis of lncRNAs in the control group (ML1-ML3) and the acute cold stress group (MC1-MC3). (B) Hierarchical clustering analysis of mRNAs in the control group (ML1-ML3) and the acute cold stress group (MC1-MC3). Yellow, high relative expression. Blue, low relative expression.

qRT-PCR validation of differentially expressed lncRNAs.

The top 10 up-regulated DE lncRNAs were selected for qRT-PCR analysis to verify the transcriptome sequencing results in all samples. The qRT-PCR results demon­strated that the 10 lncRNAs were all up-regulated in the acute cold stress group. As shown in Figure 3, the relative values of the expression levels detected by qRT-PCR were consistent with the sequencing data. This result suggests that the transcript sequencing results were reliable.

Figure 3. The expression levels of the 10 up-regulated DE lncRNAs by qRT-PCR.

Differentially expressed mRNA and lncRNAs.

The gene expression values were quantitatively estimated by using the FPKM value, and 166 lncRNAs induced by acute cold stress were screened out with the control group as the reference. Except for two known down-regulated lncRNAs, all other lncRNAs were newly discovered. A total of 269 genes showed significant changes. Among the significantly up-regulated genes, the potassium channel subfamily K member 3 (KCNK3) had the largest fold change (Log₂ Fold Change  = 7.49). Among the significantly down-regulated genes, aldolase and fructose-bisphosphate B (ALDOB) had the largest fold changes (Log₂ Fold Change = −6.24). The top 10 DE genes, including both up- and down-regulated genes, are listed in Table 2.

Table 2. The list of top 10 DE genes.

No.	Gene symbol	Description	Log2 fold change	Up/down
1	KCNK3	potassium channel subfamily K member 3	7.49	up
2	KNDC1	kinase non-catalytic C-lobe domain containing 1	5.92	up
3	TERB2	telomere repeat binding bouquet formation protein 2	5.45	up
4	CCDC38	coiled-coil domain containing 38	4.31	up
5	ZNF185	zinc finger protein 185 with LIM domain	3.99	up
6	PCK1	phosphoenolpyruvate carboxykinase 1	3.89	up
7	OASL	2'−5'-oligoadenylate synthase-like protein	3.75	up
8	CHST13	chondroitin 4-sulfotransferase 13	3.67	up
9	RHBDL3	rhomboid-related protein 1/2/3	3.44	up
10	RSAD2	radical S-adenosyl methionine domain-containing protein 2	3.27	up
11	ALDOB	aldolase, fructose-bisphosphate B	−6.24	down
12	SLIT1	slit guidance ligand 1	−5.98	down
13	MAP6	microtubule associated protein 6	−5.20	down
14	PRND	prion like protein doppel	−4.50	down
15	ANGPTL8	angiopoietin like 8	−4.42	down
16	KLHL31	kelch like family member 31	−4.12	down
17	MEI1	meiotic double-stranded break formation protein 1	−3.99	down
18	KCNH3	potassium voltage-gated channel subfamily H member 3	−3.59	down
19	LEP	leptin	−3.20	down
20	ADAMTS19	ADAM metallopeptidase with thrombospondin type 1 motif 19	−3.18	down

Construction of a lncRNA-mRNA co-expression network

A co-expression network was drawn according to the identified relationships between DE lncRNAs and mRNAs and the genes predicted by cis and trans targets of the lncRNAs. The lncRNA-mRNA network we constructed consisted of 10 up-regulated lncRNAs (blue rectangles in Figure 4) and 46 mRNAs (red or green triangles in Figure 4). Among the 10 up-regulated lncRNAs, only MSTRG.103678 had one target gene (KCNK3); all others all had between 3–19 target genes. As shown in Figure 4, MSTRG.77107 was co-expressed with 19 of the 46 mRNAs identified, the most of all identified DE lncRNAs. The next most abundant lncRNA, MSTRG.107837, was co-expressed with as many as 13 mRNAs. Moreover, many mRNAs were co-expressed with the maximum number (three) among the 10 DE lncRNAs, including BDKRB2, COL28A1, EML6, FOXN2, ITGB8, KIAA1324L, KLF6, MLPH, and PAK1.

Figure 4. lncRNA and mRNA co-expression networks in subcutaneous adipose tissues of the acute cold stress group. The co-expression network consists of the top 10 up-regulated DE lncRNAs (blue rectangles) and the correlated 46 mRNAs (red or green triangles; red represents genes that were up-regulated, and green represents genes that were down-regulated).

GO enrichment and pathway analysis

To further explore potential targets of the DE genes in the acute cold stress group, GO enrichment and KEGG pathway analyses were performed on the DE genes (Figure 5). GO analysis usually comprises three catego­ries: Biological processes (BP), cellular components (CC), and molecular functions (MF). Numerous GO terms were targeted by these coding mRNAs, including 32 biological processes (redundant GO terms removed) and one molecular function, but no cellular components (corrected P value < 0.05). For the DE mRNAs in the acute cold stress group identified by transcriptome sequencing, the most significant emerging themes as determined by corrected P-value were ‘defense response to virus’ (BP; Figure 5) and ‘double-stranded RNA binding’ (MF).

Figure 5. The enriched GO terms of ‘biological processes’ for DE mRNAs in the acute cold stress group based on the corrected P-value.

Interestingly, although there were many genes with altered expression, no pathways were enriched significantly (i.e. having a corrected P value < 0.05). The top five pathways were PPAR signaling pathway, Retinol metabolism, ECM-receptor interaction, Neomycin, kanamycin, gentamicin biosynthesis and Glycolysis / Gluconeogenesis. However, none reached the normal significance level. This also shows that acute cold stress had no significant effect on the bodies of Min pigs, consistent with the breed’s cold resistance.

Discussion

Low temperatures can change the metabolic pattern of fat deposition, especially acute cold exposure (Gao et al. 2018). In this study, the relatively cold-resistant Min pig was used as the experimental subject. The expression patterns of mRNA and lncRNA in subcutaneous fat after 24 h of acute cold exposure were analyzed by transcriptome sequencing. The number of up-regulated lncRNAs was slightly lower than that of down-regulated lncRNAs (78 vs 88), but the number of up-regulated mRNAs was slightly higher than that of down-regulated mRNAs (141 vs 128). The change trend of the two was not completely consistent.

In the analysis of differentially expressed genes, we first observed the expression changes of genes related to lipid metabolism. ANGPTL4 (Fernández-Hernando and Suárez 2020), KLF11 (Zhang et al. 2013), and NUPR1 (Teresa Borrello et al. 2021) related to lipid metabolism were significantly up-regulated by two- to threefold. However, expression levels of beige adipose tissue marker genes CD137, TMEM26, Cited1, and Tbx1 as well as key thermogenic genes UCP3, PGC1α, Dio2, CPT1α, Cidea, PPARα, and ADRB3 were not significantly different. The expression of the Leptin gene, which regulates fat storage and speeds up metabolism, was decreased by a factor of 3.2. A decrease in its expression level suggests that cold stimulates the body's food intake (Akana et al. 1999). The results indicate that cold stress did change the metabolism of subcutaneous fat in Min pigs, but there was no obvious browning.

Further analysis of other differentially expressed genes showed that multiple interferon stimulating genes (ISGS) such as RSAD2, ISG20, ISG15, DDX60, MX2, MX1, and IFI6 were significantly up-regulated. These genes play key roles in the process of host antiviral defense (Tirumurugaan et al. 2020; Zheng et al. 2017). In addition, DHX58 is related to viral infection and the number and distribution of immune cells in the blood (Li et al. 2010); PTX3 is related to innate immune response (Daigo and Hamakubo 2020); OASL is involved in resistance to classical swine fever virus (Cai et al. 2017), and HERC5 can regulate ISG15 and respond to multiple virus infections (Mathieu et al. 2021). Their transcription levels were significantly increased, indicating that cold stimulation induced the pigs’ natural immune system, a result consistent with the findings in mice (Liang et al. 2019). The results of GO annotation also suggest that Min pigs do not immediately start non-shivering thermogenesis within 24 h of cold stress, but rather initially react via the body's defense and immune response to resist cold stimulation and improve the body's immunity.

Cold stress can reduce brain blood flow (CBF) by 20–30% (Gibbons et al. 2020), resulting in poor local blood circulation, tissue ischemia and hypoxia, and even necrosis (Thomas et al. 1994). In previous studies, our team found that the ear tissue and body surface of Large White pigs become red, swollen, and necrotic after acute cold stress. Cold stress is often accompanied by increased respiration and oxidative stress (Blagojevic et al. 2011). After 24 h of acute cold stress, many genes related to blood circulation, respiration, and oxidative stress showed significant expression changes in Min pigs. For example, the KCNK3 gene had the largest fold change (log₂ fold change  = 7.49); its function loss is the key to the pathogenesis of pulmonary hypertension (Lambert et al. 2019). The CCDC38 gene is expressed in bronchial epithelial cells and is associated with one second rate (FEV1/FVC, the proportion of forced expiratory volume in the first second to all expiratory volume) (Wain et al. 2014). A variant of GAS2L2 leads to ciliary immobility syndrome that characterized by recurrent respiratory tract infection (Bustamante-Marin et al. 2019). KNDC1 (Ji et al. 2018) and LIPG (Leonhardt 2019) are related to oxidative stress; C1QL3 is related to angiogenesis, and ETV7 is related to hematopoiesis. All were significantly up-regulated (Li et al. 2019; Numata et al. 2019).

Genes related to glucose metabolism such as the rate limiting enzyme PCK1 (Xu et al. 2020) regulating gluconeogenesis or PDK4, which plays an important role in glucose consumption and fatty acid conversion (Park et al. 2018) were up-regulated, while ALDOB, which mediates fructose metabolism, was down-regulated.

Through the joint analysis of lncRNA and mRNA, this study found that only 7% of the target genes were cis regulated by lncRNA, and more target genes were trans regulated, a result that was basically consistent with the prediction results of 6.8% by Zhan et al. (2016). In the co-expression network constructed by lncRNA and mRNA, there was one mRNA that could bind multiple lncRNAs and one lncRNA that could interact with multiple mRNAs, resulting in a very complex relationship. Taking lncRNA MSTRG.90687 as an example, it could interact with six genes at the same time, UCP3, KLF11 and PLIN5 related to fat metabolism, TTC39A and ADAMTS18 related to reproduction, and SLC37A4 related to glycogen metabolism. The Leptin gene related to fat metabolism could also bind to six different lncRNAs (MSTRG.79145, MSTRG.109887, MSTRG.157809, MSTRG.27578, MSTRG.29107, and MSTRG.67063) at the same time. Because the pig lncRNA database is not complete, the functions of many lncRNAs are unknown. Which lncRNAs play an important regulatory role in the response to cold stress needs to be verified by further experiments.

It can be seen that, the mRNA and lncRNA of Min pig backfat were changed after acute cold stress. But no pathway was significant enriched, it means the physiological function of Min pig was not affected after short acute cold stress.

Availability of data and materials

All data generated or analyzed during this study are included in this published article.

Authors’ contributions

ZD, WW and LD conceived and designed the study. WW, LZ and WL collected the samples and acquired the data. ZD prepared the manuscript. All authors reviewed and approved the final manuscript.

Ethics approval and consent to participate

The study protocol was conducted in accordance with the Declaration of Helsinki and was approved by the Institutional Ethical Review Board of Heilongjiang Academy of Agricultural Sciences. All participants provided written informed consent.

Patient consent for publication

Not applicable.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The present study was Science Foundation for Excellent Youth Scholars of Heilongjiang Province (JQ2020C005)