Population Samples
In total, 186 samples of Chinese wagyu beef cattle(Both F2 and F3 cross-population of Luxi 、Qinchuan and Wagyu) were taken from three different farms in China (Shandong Kaiyuan Co. Ltd., Shandong Province, China; Ningxia Yijiayi Farming and Animal Husbandry Co. Ltd., Ningxia Hui Autonomous Region,China;Ningxia Xuanheyuan Agriculture and Animal Husbandry Co. Ltd., Ningxia Hui Autonomous Region, China). The background information of samples is shown in table 1.
Cattles are hybridized for four generations according to fig.1. The first generation is the cross between Qinchuan and Wagyu cattle, The offspring of the cow (F2) is used to continue the cross with Wagyu (All hybridization was performed by artificial insemination) for another two generations. The castration procedure was performed at 28 months of age before the cull. Blood and longissimus dorsi muscles were collected and intramuscular fat content was determined. The second and third generations were bred and raised in the same way as the first generation(fig 1).
Subject to the limitation of cattle farms hybridization program, experimental samples from one single cattle farm are very difficult to achieve. Therefore, we selected three cattle farms that were certified for wagyu importation and sperm in-vitro fertilization. The source of wagyu beef sperm came from Beijing Dairy Cow Center(BDCC). The Simmental beef cattle used as control were obtained from Ningxia Muquan Halal Slaughterhouse in China. The animal use protocol listed below has been reviewed and approved by the Animal Ethical and Welfare Committee(AEWC). NO. IACUC-NXU1014, and all cattles were electric shocked before slaughter, Electric shocks can make cattles painless when they are slaughtered, This test conform to the requirements of American Veterinary Medical Association (AVMA) Guidelines. Samples collected in this experiment have been given permission from respective farms in China from where of Chinese Wagyu Beef Cattle. I've provided three statements to confirm that permission was obtained from respective farms in China from where the of Chinese wagyu beef cattle. Three statements are uploaded in the Supplementary Materials and Related Files section.
Phenotypes determination
Phenotype determination was obtained after acid excretion and segmentation. In short, 100g of eye fillet muscle was cut out, and then intramuscular fat content was determined according to the national fat content determination standard. The classification of eye fillet meat is based on the Chinese beef grading standard, which includes the carcass and the yield quality rating. Both quality ratings were conducted after the cattle carcass was frozen and the acid was removed. The marbling of the 12th-13th rib longissimus dorsi muscles in cross-section and the top quality of the cattle and the meat and fat color was used as the reference standards. According to the amount of intramuscular fat in the cross-section, the grading of marbling is described as follows: extremely rich intermuscular fat is grade 5 and gradually reduce to grade 1 which has the least intramuscular fat. In addition to the marbling score, the samples of physiological maturity were determined by the ossification score. There is a tendency that cattle with more marbling often had less physiological maturity, which means the younger of the age and overall, the higher the beef grade score. It should be noted here that the final grade has been adjusted by the color of meat and fat and standardized across different farms.
The intramuscular fat content was determined following the guideline of Chinese national standard GB/T 9695.7-2008 in the booklet "Determination of Total Fat in Meat and Meat Products". The crude fat content was determined by Soxhlet method.
Table 1 The background information of samples
Dam
|
Male
|
Number
|
Name of the cattle farmer
|
Qinchuan
|
Wagyu
|
113
|
Ningxia Yijiayi Farming and Animal Husbandry Co. Ltd., Ningxia Hui Autonomous Region,China
|
Qinchuan×Wagyu
|
Wagyu
|
40
|
Ningxia Xuanheyuan Agriculture and Animal Husbandry Co. Ltd., Ningxia Hui Autonomous Region,China
|
Qinchuan×Wagyu×Wagyu
|
Wagyu
|
33
|
Shandong YuanLong Co. Ltd., Shandong Province, China;
|
DNA Sequencing
A total of 2 µg genomic DNA (35 ng/ml) was fragmented into 300 bp using the Bioruptor UCD-200 sonicator (Diagenode, Denville, NJ) for each sample. Library preparation was constructed using the Kapa Hyper DNA library prep kit for the Illumina platform. Fragmented DNA was end-repaired with an end-repair enzyme and a deoxyadenosine was added to the 3’ ends of the fragments. Kapa barcoded DNA and Kapa indexed adapters were ligated to the sample libraries. The adapter-ligated libraries were selected for an average insert size of 300 bp using next-generation sequencing cleanup and size selection kit (NucleoMag, Macherey-Nagel, Duren, Germany) according to the manufacturer’s protocols. The quality of libraries was assessed using the Bioanalyzer 4200 (Agilent Technologies, Santa Clara, CA). The libraries were then quantified by qRT-PCR and then sequenced were performed by Illumina Nova-seq platform to generate 150-bp paired-end reads.
Sequencing Genotypes Calling
Raw sequence reads were processed using Fastq(v0.12.4)10 to remove low-quality reads and adapters. BWA (v0.7.16)11 mem (-M -a) module was applied to align the high-quality reads to ARS-UCD1.2 cattle genome with default parameters. SAMtools (v1.9) 12 were used to sort BAM files, remove read of PCR duplications, and for sample alignment statistics. “BaseRecalibrator” and “ApplyBQSR” module of GATK (v4.0.10.1) were used for base quality score recalibration. “HaplotypeCaller” module with minimum-mapping-quality 30 for each sample GVCF. After that, raw cohort VCF was worked out with modules “CombineGVCFs” and “GenotypeGVCFs”. Bcftools13 (v1.9) was performed as variants quality filtering with "QUAL >=100". To annotate variants, ensemble gene annotations base on ARS-UCD1.2 were used for the database which was created by snpEff (v4.3T)14. The brief workflow of variants calling and annotation was shown as blue boxes in figure 2.
BSA-based Allele Frequency Analysis
The allele frequency of each population was extracted from the VCF file.
One of the challenges of sequencing-based methods to identify causal mutation sites is the high background noise resulting from a large number of potential SNPs, and sequencing/mapping errors. We developed a new sequence data analysis model, named as DFScore, which calculates the variance of allele frequency in a sliding window. When the allele frequency in a window obeys the frequency expectation of two bulk samples or the variance in each bulk sample close to or equal to zero, we treat it as the best window. The benefit of this calculation is that we can more accurately screen segments that meet genetic expectations. For instance, assumed the data is from an F2 population, the allele frequency expectation of dominant bulk should be 1/3 or 2/3 and the allele frequency expectation of recessiveness bulk should be 0 or 1. Similarly, if the data were generated with homozygous samples in the two extreme phenotype bulks, the allele frequency expectation of the two bulks will be either 0 or 1. The benefit of this calculation is that we can more accurately screen segments that meet genetic expectations.
The Power AFD can be calculated as follows:
PowerAFD = |FREQ of High IMF Group-FREQ of Low IMF Group|^4
Gene Expression Analysis
Gene expression and annotation analysis were obtained on the cattle gene expression database from Bgee (https://www.bgee.org/). In total, 112 muscle and fat tissues of cattle RNA-seq datasets were downloaded from NCBI and processed to the FASTQ format using the NCBI sratoolkit (version 2.9.6). Then low-quality reads and adapters were removed using the Fastq program (version 0.12.4,)10. Kallisto (version 0.45.0)15 was used to quantized gene expression. Ensembl genome bos_taurus.ARS-UCD1.2 and gene annotation version 100 were used as references gene. Gene expression heatmap was generated by R package “pheatmap” (version 1.0.12) with gene TPM.