The genome sequence of 24 Han Chinese patients with CD from 12 families were obtained (Table 1). Written informed consent was provided by the attending parents or legal guardians of the pediatric participants. This study was approved by the Research Ethics Committee (REC) of the Sixth Affiliated Hospital of Sun Yet Sen University (SYSU) (N0.2020ZSLYEC-006).

Using a QIAamp DNA Kit (QIAGEN, Hilden, Germany), genomic DNA was extracted from peripheral blood cells in accordance with the manufacturer’s protocol. DNA samples were quantified using a Qubit (Thermo Fisher Scientific). A total of 2 μg of each DNA sample was sent to the Beijing Genome Institute (BGI, Shenzhen, China) for WGS using the BGISEQ-500, according to the manufacturer’s guidelines. According to the manufacturer’s recommendations, the genomic DNA was briefly split by ultrasound on a Covaris E220 (Covaris) to DNA segments between 50 bp—800 bp. The fragmented DNA was then exposed to end-repair, phosphorylation, and A-tailing procedures after being further chosen to 100bp-300bp using AMpure XP Beads (Beckman Coulter, Indiana, United States). The A-tailed segments were ligated to the BGISEQ-500 platform-specific adaptors, and the ligated fragments were then purified and amplified using PCR. Finally, single-stranded DNA circles were produced by the circularization process. The libraries were sequenced using 50 bp paired-end reads on the BGISEQ-500 platform following quantification and qualifying.

SOAPnuke was used to filter the raw sequencing reads (Li et al., 2009b) (N rating >10%, low quality rating >50%, and quality rating<5) and Burrows-Wheeler Aligner (BWA v0.7.17) to align to the UCSC human reference genome (hg19) (Li and Durbin, 2009; Li and Durbin, 2010). The coordinates were sorted using Samtools (version 1.3.1) and duplicates were identified using Picard (version 1.129, http://picard.sourceforge.net) (Li et al., 2009a). Using GATK HaplotypeCaller (McKenna et al., 2010), single nucleotide substitution variants (SNV) and brief insertions and deletions (indels) were identified (McKenna et al., 2010). We used the GATK Variant Quality Score Recalibration (VQSR) that uses machine learning algorithm to filter the raw variant callset. The GATK VQSR used high-quality known variant sets as training and truth resources and built a predictive model to filter spurious variants. The SNPs and InDels marked PASS in the output VCF file were high-confident variation set. For SNPs recalibration strategy, we used the following datasets and features to train the model. (a) Training sets: HapMap V3.3, Omni2.5 M genotyping array data and high-confidence SNP sites produced by the 1000 Genomes Project. (b) Features: Coverage (DP), Quality/depth (QD), Fisher test on strand bias (FS), Odds ratio for strand bias (SOR), Mapping quality rank sum test (MQRankSum), Read position rank sum test (ReadPosRankSum), RMS mapping quality (MQ). For InDels recalibration strategy, we used the following datasets and features to train the model. (a) Training sets: Mills 1000G gold standard InDel set. (b) Features: Coverage (DP), Quality/depth (QD), Fisher test on strand bias (FS), Odds ratio for strand bias (SOR), Mapping quality rank sum test (MQRankSum), Read position rank sum test (ReadPosRankSum).

All germline SNV and indels were annotated using an in-house annotation pipeline, as described previously (Neveling et al., 2013; de Voer et al., 2013; Vissers et al., 2010). High-confidence calls (i.e., ≥10 reads, ≥5 variant reads, and ≥20% variant reads) were subsequently prioritized for variants that were non-synonymous and were absent in our in-house variant database (2,037 in-house analyzed exomes, mostly from European ancestry). Next, we removed all variants present with a MAF of >0.001 in dbSNPv138, the National Heart, Lung, and Blood Institute (NHLBI) Exome Sequencing Project database (ESP, 6503 exomes, http://evs.gs.washington.edu/EVS/), the Exome Aggregation Consortium (Lek et al., 2016) and Thousand Genome Project. Subsequently, family index patients shared non-synonymous variants that result in alterations in protein function, including protein truncation, splice site defects and missense mutations at highly conserved (phyloP ≥ 3.0) nucleotide positions, were included in our analyses. Alamut v.2.0 software (Interactive Biosoftware) and integrated mutation prediction software (align GVDV, SIFT and PolyPhen-2) (Adzhubei et al., 2010; Kumar et al., 2009; Vissers et al., 2010)packages were used for analyses of the identified variants. The prediction of splicing effects was evaluated based on five different algorithms (SpliceSiteFinder, MaxEntScan, NNSPLICE, GeneSplicer, Human Splicing Finder) through the bioinformatics tools of the Alamut v.2.0 software.

By using GWAS, we first targeted germline variations in genes at susceptibility loci known to be linked to IBD (Anderson et al., 2011; de Lange et al., 2017; Ellinghaus et al., 2016; Huang et al., 2017; Julia et al., 2014; Kenny et al., 2012; Liu et al., 2015; Parkes et al., 2007; Yamazaki et al., 2013; Yang et al., 2014). Genetic defects in innate immunity that impair intestinal bacterial sensing are linked to the development of IBD (Cananzi et al., 2021). Recent developments in molecular biology have uncovered crucial details about the genetic basis of numerous inflammatory diseases. Next to the identification of variants in known IBD GWAS genes, we searched for potential pathogenic variants in novel candidate genes using the remaining genome data of our CD family cohort. We concentrated on genes that fulfilled the following criteria while choosing these variants: (Supplementary Table S2): 1) genes with variations that caused protein truncation (such as putative frameshifts, nonsense variations, and variations at canonical splice sites), as well as non-synonymous variations with a PhyloP score of more than 3.0, were chosen; 2) the International Union of Immunological Societies Expert Committee on Primary Immunodeficiency’s collection of primary immune deficiency (PID) genes (Picard et al., 2015); and 3) genes involved in pathways implicated in IBD pathogenesis, including innate immune system, immune system, and neutrophil degranulation pathway (Belinky et al., 2015; Kanegane, 2018).

After PCR amplification, WGS-identified candidate variant of SIRPB1 was verified using Sanger sequencing. The Primer3 software program was used to build PCR primers in silico. Standard PCR procedures were used on an Applied Biosystems Dual 96-Well GeneAmp PCR System 9,700 (primer sequences available upon request). Using the software package Vector NTI, variant analyses were carried out (Invitrogen, Paisley, United Kingdom).

Candidate variant validation analysis was performed on 381 probands with IBD and 381 unrelated individuals recruited by the Sixth Affiliated Hospital of SYSU using MassARRAY (Derkach et al., 2013). All individuals recruited were peotected by the REC of the Sixth Affiliated Hospital of SYSU (N0.2020ZSLYEC-006). The phenotypes of the controls were assessed, and no known gastrointestinal or immunological findings were reported. Fisher’s exact testing was conducted, and statistical significance was set at p < 0.05 (Derkach et al., 2013).

All the microarray samples used in this study were systematically searched and downloaded from NCBI-GEO (Barrett et al., 2013) after the manual curation of the sample details. We obtained gene expression data of mucosal biopsies in CD patients and normal controls from following array data series: GSE75214, GSE36807 and GSE59071. The bioinformatics online tool GEO2R was used to analyze the mRNA expression of SIRPB1.

To study the variant type in vitro to understand its functional implications in more detail, we established THP-1 cell lines that stably expressed wild-type SIRPB1 (SIRPB1 ^wt ) and mutant SIRPB1 (SIRPB1 ^11143iG , c.1143_1144insG; p. Leu381_Leu382fs). Human SIRPB1 ^wt and SIRPB1 ^11143iG expression vectors (pEZ-M02/SIRPB1 ^wt and pEZ-M02/SIRPB1 ^11143iG , respectively) were established. According to the manufacturer’s instructions, Lipofectamine 2000 (Invitrogen; Thermo Fisher Scientific, Waltham, MA, United States) was used to transfect THP-1 cells with 0.5 µg of either pEZ-M02/SIRPB1 ^wt and pEZ-M02/SIRPB1 ^11143iG . Continuous neomycin treatment at 450 μg/mL was used to select transfectants that could consistently express the inserted vector plasmid (Thermo Fisher Scientific). The limited dilution approach was used to clone neomycin-resistant cells, which were then kept alive in medium containing neomycin.

By administering THP-1 monocytes 100 ng/mL phorbol 12-myristate 13-acetate (PMA; Sigma) for 48 h, the macrophage-like state was generated. Then, THP-1 macrophages were transfected with 0.5 µg pEZ-M02/SIRPB1 ^wt and pEZ-M02/SIRPB1 ^11143iG plasmid for 4 h. After transfection, THP-1 macrophages were incubated in complete medium for 48 h and stimulated with LPS (100 ng/mL). THP-1 cells were lysed in ice-cold lysis buffer for western blotting, and proteins were separated on 10% SDS page. The main antibodies against DAP12, p-Syk/Syk, p-Akt/Akt, p-Jak2/Jak2, and anti-GAPDH were then used to probe the membranes. According to the manufacturer’s instructions, the culture medium was suspended for ELISA assays and TNF-α, IL-1, and IL-6 ELISA kits were used to detect the substances.

Prism version 8.0 software (GraphPad) was used to perform the statistical analysis. Both the Student’s t-test and the Dunnett’s test were used to determine the significance of differences. When more than two groups were compared, a one-way ANOVA was performed. The data are shown as mean ± SE, and a significance level of 0.05 was used.

Studies of twins and familial clustering of disease clearly indicate that IBD, especially CD, is a hereditary disorder (Halfvarson et al., 2003). Therefore, we hypothesized that family-based CD is more probable to have occurred as a result of uncommon or unique variations and may play a role in the disease’s development. Thus, 24 patients with CD from 12 families in the case database of our IBD center were enrolled for WGS analysis. Table 1 provides a summary of the clinical traits of the individuals that were enrolled. Among them, there were two pairs of fraternal twins, three pairs of identical twins, four pairs of father or mother and daughter or son, two pairs of elder sister or brother and younger brother, and one pair of aunt and niece (Figure 1). All the patients were of Han origin, including 13 men (54.17%) and 11 women (45.83%). The median age of participants in this study was 22 years (range 14–45 years). Regarding disease location, 20 patients were diagnosed with ileocolonic inflammation (L3, 83.33%), one patient with colonic inflammation (L2, 4.17%), 3 patients with ileal inflammation (L1, 12.50%), and six patients had concomitant upper gastrointestinal disease (25.00%). Regarding disease behavior, nine patients were accompanied with non-stricturing non-penetrating lesions (B1, 37.50%), six patients had stricturing lesions (B2, 25.00%), and nine patients had penetrating lesions (B3, 37.50%). A total of 13 patients (54.17%) had pelvic disease and six patients (25.00%) had a history of abdominal surgery.

To further understand the genetic basis of IBD, we used WGS to look for uncommon, possibly disease-causing coding variations in familiar CD patients (Table1). A substantial percentage of coding variants discovered across individuals in a family were prioritized for further exploration using a series of progressive variant filters. We identified on average 3,430,329 variants (range: 3,393,972–3,466,204) per genome. A prioritization scheme was applied to identify candidate variants as shown in table 2. After rigorous bioinformatics analysis steps (see Methods), there were 92 genetic variants significantly associated with CD in Chinese individuals including several known IBD related gene including NOD2, ZPBP2, TNFSF15, LSP1, SDCCAG3, MUC19. Then, we adopted the MassArray Analyzer system (Sequenom, Inc., SanDiego, CA, United States) for genotyping those genetic variants. We designed primers for 89 variants in Sequenom official websites, while, there is suitable primer for the remaining 3 variants. However, those 89 primers were distributed in 5 panels, due to the limited budget, only the top 2 largest panels (a total of 61 variants were included, about 30 variants per panel) used for the subsequently assay. We finally generated a candidate variants list including 61 variants for replication analyses with an independent cohort including 381 patients with CD and 381 control subjects (Table 3).

Notably, we identified the frameshift variant SIRPB1 p. Leu381_Leu382fs in two index patients, F7-II-1 and F7-II-2, from Family 7 (Figure 2). Patient F7-II-1 was the elder sister of patient F7-II-2, whose age was 23 years with a disease duration of about 2 years, and patient F7-II-2 was the younger brother whose age was 21 years with a disease duration of approximately 1 year. Both manifested as ileocolonic inflammation (L3). The disease behavior of F7-II-1 was non-stricturing non-penetrating (B1), whereas the disease behavior of F7-II-2 was penetrating (B3). For replication analysis, the blood sample of 384 CD patients including 104 females and 280 males in our database were used for microarray analysis to verify the susceptibility gene mutation. The mean age of these patients was 27.0 ± 10.4 years-old. The 384 blood samples of control group were collected from healthy testing population including 170 female and 214 female and the mean age of these patients was 34.7 ± 9.1 years-old. The unique uncommon frameshift variant in SIRPB1 was considerably enriched in CD patients compared to controls, according to genotyping analysis (p = 0.03, OR 4.59, 95% CI 0.98–21.36, 81.82% vs. 49.53%). (Table 3). The clinical characteristics of patients with CD identified to harbor the SIRPB1 p. Leu381_Leu382fs variant are listed in Table 4.

A member of the family of signal-regulating proteins (SIRP) and of the immunoglobulin superfamily, SIRPβ (also known as CD172b) is encoded by SIRPB1 (van den Berg et al., 2005). Previous research on SIRPB1 mostly focused on its biochemical properties and functions and discovered that it stimulates DAP12 and Syk tyrosine phosphorylation, which then activates the mitogen-activated protein kinase (MAPK) pathway to increase phagocytosis in macrophages (Hayashi et al., 2004). However, the role of SIRPB1 in IBD pathogenesis has not yet been reported. Therefore, we first assessed the expression levels of SIRPB1 by analyzing publicly available data from patients with IBD and healthy individuals. According to the findings, patients with active CD (A-CD) had considerably higher levels of SIRPB1 expression in their ileocolonic tissue than either healthy individuals or patients with CD that was in remission (R-CD) (Figures 3G–I), indicating that SIRPB1 is possibly involved in the progression of CD. Next, we assessed the expression levels of SIRPβ in the ileocolonic tissue of variant and wild-type patients with CD using IHC (Supplementary Table S3). Notably, patients with CD and the SIRPB1 variant exhibited significantly higher levels of SIRPβ than those with CD and wild-type SIRPB1 (Figures 3A,D). Consistently, SIRPβ-mediated transduction signal molecules, such as DAP12 and p-Syk, were also significantly upregulated in patients with CD and the variant compared with WT controls (Figures 3B,C,E,F). Our data suggest that the higher expression of SIRPβ was caused by the frameshift variation of SIRPB1.

As documented above, the frameshift variant in SIRPB1 led to higher expression of SIRPβ; however, the function of the rare SIRPB1 frameshift variant still needs to be elucidated. The interaction of SIRPβ with the activating adaptor protein DAP12, which carries an immunoreceptor tyrosine-based activation motif (ITAM) and transmits activating signals, depends on the presence of a basic amino acid side chain in the transmembrane domain of SIRPβ (Lanier and Bakker, 2000; Tomasello and Vivier, 2005). An earlier study showed that phosphorylating Syk and MAPK and crosslinking mouse SIRPβ with monoclonal antibodies increases neutrophil migration and macrophage phagocytosis (Hayashi et al., 2004; Liu et al., 2005). As shown in Figure 4, the frameshift variant of SIRPB1 gene led to an alteration of the amino acid sequence located both on the transmembrane and the cytoplasm. Therefore, we hypothesized that the frameshift variant of SIRPB1 could make functional contributions to the inflammation process.

To investigate the potential effects of the identified novel frameshift variant of SIRPB1, Plasmids carrying either SIRPB1 ^wt and SIRPB1 ^11143iG sequence were transfected into THP-1 cells that had been stimulated with PMA for 48 h. In line with a previous report, compared with THP-1 cells transfected with the empty vector (NC), THP-1 cells with SIRPB1 overexpression (SIRPB1 ^wt ) exhibited increased activation of DAP12 and NF-κB and elicited tyrosine phosphorylation of Syk, Akt, and Jak2, causing increased IL-1β, TNF-α, and IL-6 release (Figures 5A,B,C,D,E). Furthermore, compared with SIRPB1 ^wt cells, THP-1 cells expressing variant SIRPB1 (SIRPB1 ^11143iG ) displayed higher SIRPβ expression, enhanced activation of DAP12 and NF-κB, increased tyrosine phosphorylation of Syk, AKT, and Jak2, and higher secretion of IL-1β, TNF-α, and IL-6 (Figures 5A,B,C,D,E). Taken together, our results show that the frameshift variant of SIRPB1 amplified the function of the wild-type SIRPB1 gene.