Total 8,295 cirrhotic portal hypertensive patients, diagnosed by medical history, signs and imaging of portal hypertension according to chinese guidelines on the management of liver cirrhosis (Xu et al., 2020), were enrolled in cohort 1. The etiologies of liver cirrhosis included hepatitis B (65.04%), hepatitis C (15.79%), alcoholic cirrhosis (12.35%), and autoimmune liver disease (6.82%). Two hundred sixty-seven patients with NCPH in cohort 2, confirmed by radiologists, hepatologists and pathologists based on enhanced computed tomography (CT) or nuclear magnetic resonance imaging (MRI) and/or pathology according to guidelines on the management of NCPH of EASL (Khanna and Sarin, 2014). The signs and symptoms, laboratory and endoscopy data were obtained from Electronic Medical Record Management System (EMRS). The portal and splenic vein diameter and splenic thickness were measured using abdominal ultrasonography. Model for end-stage liver disease (MELD) and Child-Pugh score were assessed for the severity of cirrhosis. Transient elastography liver stiffness measurement (LSM) was performed using FibroScan™ (Echosens, Paris, France). The MRC were defined as more than two cysts in both kidneys detected by ultrasonography or CT.

The study protocol was performed in compliance with the Declaration of Helsinki and approved by the Ethics Committees of Beijing You’an Hospital Capital Medical University. Signed informed consent was obtained from each participant for using samples, materials and publication.

The peripheral blood of nine patients with MRC (four from cohort 1 and five from cohort 2 separately) were collected for WGS to explore potential genetic modifiers or candidate disease genes. First, the genomic DNA was extracted from the peripheral blood of those patients. The whole-genome DNA Sequencing libraries were prepared according to the manufacturer’s instructions. The raw reads were produced by a BGISEQ-500 sequencer at an average depth of 40×. The rare putatively pathogenic variants were validated by Sanger sequencing.

Sequencing data were quality controlled with adapter and aligned to human reference genome build hg19 (http://www.gencodegenes.org/releases/19.html) with BWA aligner (Li and Durbin, 2009). The GATK best practice Haplotype Caller pipeline was implemented for SNV and indel calling (Li et al., 2009; McKenna et al., 2010). SV was called with lumpy software, and CNV was detected with FreeC software (Layer et al., 2014). All SNV variants were annotated using ANNOVAR for bioinformatics analysis (Wang et al., 2010). Several genomic databases, including the 1,000 Genomes (1000G), ExAC (Exome Aggregation Consortium), Exome Sequencing Project (ESP), gnom AD (both WES and WGS databases), and CG46, were used to assess the variant frequency in the population. MCAP, SIFT, Polyphen2-HDIV, Polyphen2-HVAR, MutationTaster, MutationAssessor, and Clinvar were implemented to annotate the effect of missense variants. GERPs were used to evaluate the conservation of the variant locus. Rare putative pathogenic variants were filtered as follows (Gao et al., 2020): the allele frequency of the candidate variant should be lower than 0.01 among 1000G, ExAC, ESP 6500, Genome Aggregation (GA) and Complete Genomics 46 (CG46) Databases (Khanna and Sarin, 2019); amino acid changing variants were kept, and GERP scores should be higher than 2.0 (Vilarinho et al., 2016); truncating variants were kept, and for missense variants, MCAP scores higher than 0.6 were automatically kept, whereas for other variants, the effect should be annotated as “Deleterious” or “Highly pathogenic” by at least two software programs for MCAP scores between 0.025 and 0.6. PKHD1 protein domain prediction was obtained using SMART (Letunic et al., 2021). Cilia genes were obtained from Syscilia.org for downstream analysis.

Five sets of scRNA-seq data were used in this study. Summarized gene expression matrices derived from multiple organs from human fetuses were obtained from the GEO database via accession number GSE156793 (Cao et al., 2020). Single-cell expression data from kidneys were generated from mice under accession number GSE140023 (Conway et al., 2020). The adult human kidney and human liver single cell transcriptome was achieved by accession GSE114530 (Hochane et al., 2019), GSE131685 (He et al., 2020) and GSE159929 (Liao et al., 2020), respectively. The single cell RNA sequencing analysis was implemented with Seurat package (Version 3.9.9). Cells were discarded according to the following criteria (Gao et al., 2020): cells that had fewer than 400 genes (UMI >0) (Khanna and Sarin, 2019); cells that had fewer than 600 UMI or over 10,000 UMI; and (Vilarinho et al., 2016) cells that had more than 15% mitochondrial UMI counts. After the above quality control, for mouse kidney, human fetal kidney, human adult kidney and human liver scRNA-seq analysis, we performed log-normalization with the “vst” method and identified 2000, 3,000, 3,000 and 3,000 variable features, respectively. We then scaled by setting the parameter “vars.to.regress” to “percent.mito” and “nCount_RNA”. Principal component analysis (PCA) was performed using the “RunPCA” function. The number of PCs was chosen by visualization plot with the “ElbowPlot” function. A shared nearest neighbor (SNN) graph was constructed using the “FindNeighbors” function with the top 40 PCs, then cells were clustered by the “FindClusters” function with the “resolution” parameter set to 0.5. The “RunUMAP” function was used for the visualization plot with the “umap-learn” method, setting “n.neighbors” to 40L, “dims” to 1:40, and “min.dist” to 0.3. Marker genes for each cluster were detected using the “FindAllMarkers” function, setting the parameter “min.pct” to 0.3 and “logfc.threshold” to 0.6. Subsequently, cell clusters were annotated manually to the major cell types according to known markers. Any cluster with multiple markers of two types of cells was manually discarded as a doublet.

The full-length human CRB3 was not annotated in the PDB Database; thus, the intact wild-type and mutated CRB3 protein sequences were annotated by Phyre2 (Kelley et al., 2015). Then, protein conformational alteration was predicted using Chimera software (Pettersen et al., 2004).

Microsoft R Open (version 3.6.1, https://mran.microsoft.com/) was used. The R packages ggplot 2 (version 3.1.0) and pheatmap (version 1.0.12) were used to generate graphs of the data. Continuous variables were compared using independent T test if data were normally distributed or Mann Whitney U test. The categorical variables were compared using X ² tests performed with SPSS 22.0 (IBM, United States).

The prevalence of MRC in NCPH patients accounted for 19.5% (52/267). It was significantly higher than that in cirrhotic patients with portal hypertension (6.2%, 513/8,295), p < 0.05. The clinical characteristics of NCPH patients with MRC were compared with 92 hepatitis B virous (HBV) related cirrhosis without MRC, which randomly selected (1:2) in cohort 1 (Table 1). The NCPH patients with MRC had a relatively older onset complications of portal hypertension (p < 0.05). In terms of the manifestations of portal hypertension, the proportion of EGVB was prominent (p < 0.05), while ascites was less (p < 0.001). The platelet counts were higher than that of hepatitis B cirrhosis (p < 0.001). Although there was no significant difference in LSM between the two groups, the ratio of platelet counts/LSM in NCPH patients with MRC was significantly higher than it in hepatitis B cirrhosis without MRC patients (p < 0.001). The clinical characteristics of nine patients with MRC enrolled for WGS were summarized in Table 2. Three patients (P1-3) showed CHF by liver biopsy (Figure 1A), in which two patients exhibit early-onset in childhood or adolescence and underwent splenectomy procedure before their visit to our hospital. The endoscopy and CT imaging confirmed EGV, dilated bile ducts and multiple renal cysts in those patients (Figures 1B,C).

The others NCPH patients (P4 and P5) showed late-onset of complication, their clinical manifestations, including large hepatorenal cysts and normal LSM, implied underlying HFDs although there is no pathological evidence. Moreover, four HBV-positive patients (P6-P9) were diagnosed at middle age with normal liver function may also pinpoint to the need of dissecting genetic factors for HFD phenotypic expression.

The copy number variations or structural variants spanning known HFD-related genes were not observed. Subsequently, we identified all gene mutations in the known Caroli syndrome-causing gene PKHD1 in 3 CHF patients with NCPH (Supplementary Table S1). All missense mutations in PKHD1 had GERP scores higher than 5.4, indicating that these mutations were in evolutionarily highly conserved regions. Moreover, all three patients harbored an additional truncation mutation, which were not reported in any public databases. Therefore, we assumed that the patients carried recessive mutations in PKHD1, and validated the mutations by Sanger sequencing (Figure 2A). In addition, analysis of PKHD1 mutation distributions in different protein domains showed that all mutations were located in the extracellular domain (Figure 2B). Further, by using scRNA-seq data to explore PKHD1 expression, we found that PKHD1 was highly expressed in ureteric bud cells and stromal cells and moderately expressed in metanephric cells in the kidney. In contrast, PKHD1 was highly expressed on hepatoblasts in the liver (Figure 2C). Moreover, we also detected PKD1 compound mutations, which may explain that the patient P1 had early disease onset at the age of five and large renal cysts (38 mm).

To narrow the potential HFD phenotype-associated genes, we retrieved all known cilia genes from the literature and the European project SYSCILIA (Boldt et al., 2016). After prioritizing candidate disease-associated cilia genes by scRNA-seq analysis, the cilia-related mutations derived from the enrolled patients were listed in Supplementary Table S1, such as CRB3 p.P114L. We first plotted the expression for candidate cilia genes in multiple human fetal organs, and annotated different cell types, particularly in the kidneys and liver (Figures 3A,B). We further assessed the cilia genes of cell clustering and annotated different cell type in adult mice kidney and adult human liver using uniform manifold approximation and projection (UMAP). We found that CRB3, TUBA4A, PTCH1 and CEP290 co-expressed with PKHD1 (Figures 3C,D). Similar co-expression pattern were also seen in human fetal kidney (Figures 3E,F) and adult kidney (Figures 3G,H). Followed by cell clustering and annotation in adult human liver, the CRB3 and PKHD1 co-expression were spotted in hepatoblasts (Figures 3I,J). To further investigate the functional effect of the CRB3 p.P114L variants, we performed protein structure remodeling and predicted that the mutation would lead to protein conformational alterations in the PDZ-domain (Figure 4A). The bona fide tight junction assembly acquired the capacity for PRKC1/PARD6A complex translocation to the apical surface by interacting with the CRB3 C-terminus. Hence, we speculated that this mutation might ultimately blocking CRB3 function. With additional scRNA-seq data, we believe that this newly proposed analytical strategy may help clinicians and geneticists to map disease-related genes (Figure 4B).