Biomarkers and Applications (ISSN: 2576-9588)

Article / research article

"Characterizing a Focused Landscape of Familial Acute Respiratory Distress Syndrome"

Inimary T. Toby1*, Neal J. Thomas2, Nithyananda Thorenoor2, Debbie Spear2, Susan DiAngelo2, Joanna Floros2,3*

1Department of Biology, University of Dallas, Irving, Texas, USA

2Center for Host defense, Inflammation, and Lung Disease (CHILD) Research, Department of Pediatrics, USA

3Obstetrics & Gynecology, The Pennsylvania State University College of Medicine, Hershey, Pennsylvania, USA

*Corresponding author(s): Inimary Toby, Department of Biology, University of Dallas, Irving, Texas, USA
Joanna Floros, Center for Host defense, Inflammation, and Lung Disease (CHILD) Research, Department of Pediatrics and Obstetrics & Gynecology, The Pennsylvania State University College of Medicine, Hershey, Pennsylvania, USA

Received Date: 07 September, 2020; Accepted Date: 11 September, 2020; Published Date: 18 September, 2020


Acute Respiratory Distress Syndrome (ARDS) affects approximately 190,600 patients per year in the United States, with mortality up to 45%. ARDS can occur as primary disease due to various factors (e.g. bacterial or viral pneumonia, lung contusion and toxic inhalation) or as secondary disease due to sepsis, pancreatitis and severe trauma. We hypothesized that ARDS-affected individuals have patterns of variants in their physiological repertoire that can be tracked, identified and utilized clinically. The goals of this study were to: (1) characterize the landscape of variants within protein coding and UTR regions in ARDS using an Exome sequencing approach; and (2) determine the variations in signaling pathways and identify functional consequences of ARDS. Towards this, we assessed an ARDS-affected individual, GP7, in the context of unaffected family members as well as unrelated individuals with ARDS, in order to elucidate underlying inheritance patterns of “private variants”. Private variants consisted of variants shared by ARDS cases but not found in unaffected individuals. Whole exome sequencing yielded 3,516 variants (represented by 2,354 genes), that were highly enriched in the primary case, GP7. Of these, there were 128 private variants and ~19% of these represented novel variants derived from this study. Variants identified and subsequent gene mapping analysis, demonstrate that there are important biological pathways and functions that distinguish ARDS cases from non-ARDS individuals. These include cell death and survival and liver hyperplasia/hyper proliferation. These in-silico discoveries highlight novel variants shared by ARDS cases, which can be further explored in deeper learning of ARDS.


Acute respiratory distress syndrome; Biological pathways; Exome sequencing; Variants


Acute Respiratory Distress Syndrome (ARDS) is a syndrome of hypoxic respiratory failure characterized by diffuse pulmonary infiltrates and accumulation of protein-rich pulmonary edema that cause reduction in lung compliance alveolar collapse and ventilation-perfusion mismatch [1-6]. ARDS affects approximately 190,600 patients per year in the United States, with mortality up to 45% [7]. Despite improvements in intensive care during the last fifteen years, ARDS is still the major cause of mortality and morbidity in intensive care [1,2,5-7]. In fact, ARDS therapy has seen limited progress since its initial description in 1967 and management is still largely supportive, with no established therapies targeted at the primary disease processes [8]. Accordingly, there is a need for methods of early detection [9]. There has been recent recognition of the clinical and biological heterogeneity within ARDS [10-12] that reflects our incomplete ARDS 3 3 understanding of the biology of ARDS. Additional contributions to the knowledge about inheritance of ARDS and/ or pathogenesis will be of great benefit in moving forward with successful clinical translation of new diagnostic, preventive, and therapeutic strategies [13-16].

ARDS occurs within one week of a known clinical insult or after worsening of respiratory symptoms. It is a consequence of various risk factors including direct (e.g., bacterial or viral pneumonia, gastric aspiration, lung contusion, toxic inhalation, and near drowning) or indirect (e.g., sepsis, pancreatitis, severe trauma, massive blood transfusion, and burn) lung injury [1-6,17-19]. There is little knowledge on the temporal relationship between detectable inflammatory changes and the onset of increased lung density, which is the current radiographic diagnostic marker for ARDS. A better understanding of these key temporal and topographic processes may contribute to advance diagnostic and prognostic biomarkers and effective therapies [7,20-22]. Genome sequencing studies on ARDS have concentrated on identification of plasma biomarkers that may facilitate diagnosis of ARDS which could, in theory, improve clinical care, enhance our understanding of pathophysiology, and be used to enroll more homogeneous groups of patients in clinical trials of new therapies, increasing the likelihood of detecting a treatment effect [2,23-28].

Most recently, exome sequencing studies have been reported for ARDS as part of an outgrowth of the NHLBI’s Exome Sequencing Project [29,30]. A potential limitation of these studies is that they have rarely been conducted on families with sample collection across several generations. An advantage of exome sequencing is that it allows for the analysis of “private” gene variants-variants that may have arisen de novo in one individual or family and thus would not be detected in another [30]. We hypothesized that ARDS-affected individuals have patterns of variants expressed in their physiological repertoire that can be tracked and utilized to complement approaches to clinical diagnosis and/or clinical monitoring. To address this hypothesis, we utilized an exome sequencing approach, which focuses on just the protein coding but also UTRs sequences in a given sample. Our data indicates that, there are unique variants and signaling pathways in ARDS cases which differ from those observed in unaffected individuals; and that the variant expression patterns observed in the familial cohort are markedly different from that of unrelated ARDS cases.

Materials and Methods

Study information and Sample Collection

Family members and control cohort were recruited, enrolled after informed consent under a protocol approved by the Human Subjects Protection Office and Institutional Review Board from the Pennsylvania State University College of Medicine. Three sets of samples were collected: primary ARDS case, and related family members; a third set of samples from unrelated ARDS subjects used in the present study were collected as previously described [31,32]. Together, we analyzed a total of 18 samples.

Whole Exome sequencing analysis and data collection

The workflow in Figure 1 outlines the key steps involved in sample data collection, data filtering and the downstream data analysis processes that we applied to the sequencing output in order to perform comparisons. Variants for each sample were identified based on the GATK [33] best practices pipeline. The bwa v0.7.3a software were used to align the paired end exome sequences to the hg19 reference and the Picard v1.102 Mark Duplicates tool was used to remove duplicates. Local realignment around indels was performed by running the GATK Realigner Target Creator and Indel Realigner tools, using the Mills and 1000G Gold Standard indels as the known indels. Base quality score recalibration was performed using the GATK Base Recalibrator tool with dbSNP build 138 and the Mills and 1000G Gold Standard indels as known sites, followed by GATK Print Reads. Variants were called using the GATK Haplotype Caller tools with the following parameters: ERC GVCF, LINEAR variant index type and 128000 variant index parameter followed by the GATK Joint Genotyping tool. The ANNOVAR v2015-03-22 [34] was used to functionally annotate the genetic variants, including when applicable gene membership (e.g. intron), conservation, nonsynonymous amino acid substitution, SIFT prediction, Polyphen2 prediction, etc.

Filtering and Classification of Variants

Ingenuity Pathways Analysis (IPA) was used to determine genes associated with ARDS from the curated literature, and to identify significantly enriched canonical pathways, networks, diseases and biological functions and upstream regulators from amongst filtered list of variants [35]. For our first step in identifying enriched variants, we extracted variants based on the primary case (GP7), control samples from the family members (GP3-5 and GP8) and ten unrelated ARDS cases, (JF9-18). While samples GP1, GP2, and GP6 were disease free, we excluded them from our controls since these originated from younger individuals and we were focused on late onset of ARDS during adulthood. However, these samples were later utilized in all follow-up clustering analysis performed to assess variant inheritance patterns. Enriched variants were filtered using the following set of criterias: (1) that GP7, which was the primary case in the family samples, had the variant; (2) at least 2 of the reduced (n=4) family controls didn’t have the variant; (3) at least 2 of the ten (JF 9-18) samples had the variant and; (4) that there were at least twice as many unrelated cases with the variant than without the variant. Using this criteria, we generated the list of variants enriched in our case sample, GP7, which were represented across 16 gene/exonic functional features. Next, these enriched variants were then classified into functional subsets for further analysis. The variants were categorized according to their respective genes and placed into separate groups. They were grouped from A-D, based upon the following criteria: (1) Group A comprised of literature genes curated from ARDS-related studies currently available in IPA; (2) Group B comprised of variants represented by exonic genes affecting stop codon, genes with frameshift mutations and non-synonymous mutations; (3) Group C- comprised of variants represented by genes (which included all genes from Group B) plus genes found within the 3’ UTR, 5’ UTR, and non-coding RNAs and; (4) Group D consisted of variants that were represented by genes shared by GP7 and all JF unrelated ARDS cases from this study. IPA analysis was performed for each group of genes.

Statistical metrics applied to pathway analysis

IPA [35] applies statistical assessments to determine pathway relationships for a given list of genes. To apply the underlying statistical metrics for a given dataset, we first uploaded a list of genes and performed a Core Analysis with the default settings in IPA. The Canonical Pathway Analysis in IPA associates the genes with the canonical pathways in Ingenuity’s Knowledge Base and returns two measures of association: (1) a ratio of the number of genes from the list that maps to the pathway divided by the total number of genes that map to the same pathway and; (2) a p-value of the Fisher’s exact test. To identify ARDS-specific significant gene sets that were within the top scoring targets from pathway assessments of the list of genes, we examined the top pathways from the following five IPA analysis modules: (1) canonical pathways; (2) upstream regulators; (3) diseases and bio functions; (4) Toxicity (tox) functions and; (5) networks. A Bonferroni correction (i.e. p value < 0.05/25 = 0.002) was applied to the p-values in order to maintain the type I error at 5%. For the top gene sets obtained from the query gene list, we kept genes with consistent empirical and biological relationships. The biological relationships between genes in a canonical pathway are referred to as IPA pathways. Using the biological relationships from the IPA pathway and the top genes in the canonical pathway as references, we derived a correlation trend between ‘Within Patient Expression Changes’ (WPEC) and ‘ordered categorical Multiple Organ Failure’ (ocMOF) labels that were biologically driven for all genes in the same pathway [23]. If this trend is consistent with the one computed from the data, the gene is retained and used to compute the dominant trajectory. The IPA pathway provides a graphical representation of the biological relationships between genes in a canonical pathway, where nodes represent genes and edges represent the biological relationships. Each edge is supported by at least one reference from the literature, a textbook, or canonical information stored in the Ingenuity Pathways Knowledge Base, providing us a relationship summary between genes. A -log p-value of 1.30 represented the cutoff for statistical significance (i.e p-value<0.05) thus any pathway with a value less than this, was filtered out of our final dataset.

Clustering of variants and heatmap assessments

The filtered variants were extracted and parsed for all samples from the study. These tables were then imported into an excel spreadsheet and calculations were done for the presence (assigned a 1) or absence (assigned a 0) of variants within a specific sample using VB scripting [36]. The final csv file was uploaded into R [37]. We then applied the ‘Heatmap’ and ‘clustergram’ scripts using R packages g plots g data, g tools, and rcolorbrewer to render a 2-d color image of the data showing the samples on the x-axis. To organize these data and identify potential relationships among presence/absence of specific variants, we utilized a hierarchical clustering with Euclidean distance metric and average linkage to generate the hierarchical tree. This type of clustering enabled us to find the similarity or dissimilarity between every pair of objects in the data set, group the objects into a binary, hierarchical cluster tree, and determine where to differentiate the hierarchical tree into clusters.

Clustering of genes and annotation analysis

DAVID (Database for Annotation, Visualization, and Integrated Discovery) is a Web-based application that provides a high-throughput and integrative gene functional annotation environment to systematically extract biological themes behind large gene lists. High-throughput gene functional analysis with DAVID helps to provide important insights that allow investigators to understand the biological themes within their given genomic study [38-40]. We took gene lists from groups: (B), (C) and (D) (Table 2) and performed clustering analysis using the DAVID tool. We were unable to perform this type of analysis for group A, as the curated gene list was proprietary information within the IPA software and is retained by Agilent Biotechnologies. We performed all searches using the default parameters, as referenced in the software manual, based upon the “highest” classification stringency cutoff to obtain functionally related gene groups. For the annotation analysis generated from DAVID, we searched all human sequence related metadata available in the database and we then extracted the data to explore gene cluster relationships for each of our gene lists. DAVID uses a set of fuzzy classification algorithms to group genes based on their co-occurrences in annotation terms and ranks the gene groups using an internal (EASE) score [41].


Variants identified from a family with one ARDS case and unrelated ARDS cases

Figure 1 shows the workflow for data collection, data filtering, and data analysis utilized in the study. The family pedigree is shown in Figure 2. The individual that was the primary focus of the study is denoted as GP7. GP1, GP2, and GP6 were the youngest family members. Variants from each sample were identified and their frequency was calculated the primary case, GP7, comprised of 3,516 highly enriched variants that were represented by 2,354 genes. A summary of the annotated functions from the list of variants yielded the following characteristics: 343 variants were exonic non-synonymous, 8 exonic frameshift, 16 exonic non frame shift, 356 exonic synonymous, 16 exonic unknown, 66 ncRNA exonic, 2 exonic stop gain, 118 in the 5’UTR region, 782 in the 3’UTR, 1343 were intronic/ncRNA intronic, and the remaining variants were categorized as having upstream, downstream or intergenic functions (Table 1).

Heatmaps were generated from a ranked frequency occurrence assessment from the list of variants found to be shared between GP7 and each of the GP family samples (Figure 3A) or JF unrelated ARDS samples (Figure 3B). Figure 3A illustrates the variant expression patterns and clustered relationships between GP7 and each of the GP family members. Whereas, Figure 3B illustrates the variant expression patterns between GP7 and each of the ten unrelated ARDS cases (JF9-18). The red color in each of the heat maps denotes a presence of the same variant as found in GP7, and the green color denotes absence of the variant. Compared to the GP7 the next highest frequency of variants found amongst the family cohort was in the GP1 sample, which represented one of the younger family members (Figure 3A). This was accompanied by hierarchical clustering, which identified GP1 as the closest in variant expression pattern to GP7 (Figure 3A).

Additional hierarchical clustering assessments showed that there were 2 clades within the family cohort. The first clade was comprised of GP7, GP1, and GP4. The second clade was comprised of GP3, GP6, GP8, GP5 and GP2. The clustering profiles indicate that the most similar sample to GP7 from the second clade was GP2 (Figure 3A). Furthermore, compared to GP7, JF14 was the highest scoring sample in terms of frequency of variants present as compared to any of the other unrelated controls (Table 2). Quantification of the variant occurrence frequency showed that JF14 was ~85.5% similar to GP7 (Table 2). Clustering profiles also demonstrated that JF14 was the most similar to GP7 in variant expression pattern, thereby it was clustered closest to GP7 (Figure 3B). For the unrelated ARDS samples, JF10 had the lowest detected frequency of variants and was clustered to the farthest left of the heatmap with respect to GP7 (Table 2 and Figure 3B). The variant expression patterns show that JF10, JF16, and JF18 were the least similar to GP7, as depicted by the green colored regions denoting absence of the variant. As a result, these were all clustered to the farthest left (Figure 3B). Interestingly, heatmap assessments that included both family members and unrelated ARDS cases, showed that GP7 is most similar, in variant expression, to the unrelated ARDS cases as compared to the family members with the exception of JF10 (S1 Figure). Further comparison of the GP7 case and the unrelated ARDS cases showed that there were 128 variants shared between GP7 and all unrelated ARDS, which were represented by 104 genes (Table 3). From this 104 genes, we found that ~19% of these were unique variants comprised of 24 variants represented by 9 genes (Table 4) and were not found in any of the other groups shown in Table 3. 3 of the 128 variants were also found in group B and were represented by the MYH14 gene.

Prediction of ARDS-related biological pathways and functions from ARDS case enriched variants

We extracted all pathway predictions for gene lists from groups (A), (B), (C) and (D) (Table 3). These were then compiled to review the statistical outcomes and determine unique pathways (S1 File). Variants were assessed for presence or absence in group (D) vs (A), (C) vs (A), and (B) vs (A), (S2 File). The outcomes from each assessment was a list comprised of variants identified along with their function and exonic information where appropriate (Additional file 2). Next, using the filtered genes from the GP7 case, we observed a ~2-fold and ~4.8-fold increase in the pool of ARDS genes for group (B) vs (A) and (C) vs (A), respectively (Table 3, Figure 4). Further assessment for each group individually based upon the presence or absence of statistically significant pathways and functions revealed one novel function, liver hyperplasia/hyper proliferation. The total number of genes found to be implicated in this function from (B), (C), and (D) was 387, 139, and 63 respectively (Figure 5). To further place this finding in context of the overall number of genes present from each gene list, ~61% of the genes from list (D), as compared to 44% from (B) and 39% from (C) were implicated in this function (Figure 5, S1 File). This function was not detected in group (A) which was comprised of genes found in the ARDS literature (S1 File).

Pairwise assessment of the biological pathways for genes derived from groups (B), (C), and (D) as compared to (A), showed that there were 6 significant pathways for (B), namely: Cellular Effects of Sildenafil (Viagra), CCR3 Signaling in Eosinophils, Inhibition of Matrix Metalloproteases, Actin Cytoskeleton Signaling, Choline Degradation I, and PRPP Biosynthesis I and none of these were previously identified as significant pathways in (A) (S3 File). There were 12 significant pathways for (C) and 5 of these (Hepatic Fibrosis / Hepatic Stellate Cell Activation, Acute Phase Response Signaling, Death Receptor Signaling, Osteoarthritis Pathway, Airway Pathology in Chronic Obstructive Pulmonary Disease) were also predicted as statistically significant in (A) (S1 and S3 Files). The other 7 significant pathways (Choline Degradation I, RhoA Signaling, ATM Signaling, Inhibition of Matrix Metalloproteases, Role of PKR in Interferon Induction and Antiviral Response, Oxidative Ethanol Degradation III, and Glutamine Biosynthesis I), were not previously identified as significant pathways in (A). There were 7 statistically significant pathways for (D) and only 1 of these (Tight junction signaling) was also predicted as statistically significant in (A) (S1 and S2 Files), whereas the other 6 (Acyl-CoA Hydrolysis, Ceramide Biosynthesis, Actin Cytoskeleton Signaling, PTEN Signaling, Formaldehyde Oxidation II (Glutathione dependent), and Epithelial Adherens Junction Signaling) were not. There were no pathways shared between (B), (C) and (D) as compared to A. However, there were 3 pathways shared by (B) and (C) (Inhibition of Matrix Metalloproteases, Actin Cytoskeleton Signaling, and Choline Degradation I). Three-way pathway interaction assessments resulted in the discovery of 1 biological process that was shared by group (A), (B), and (C) but was not present in group (D). This biological process contained genes with functions in cell death and survival. We found there were 95, 123 and 237 genes implicated in this function from group (A), (B), and (C), respectively (Figure 6). Of the 95 genes from (A), 82 were shared by (A), (B), and (C) (S3 File). The number of additional genes contributed from (B) and (C) was 41 and 155, respectively, and this correlated positively with the increase in the gene search pool (Figure 6, S4 File). Additionally, we found 3 variants (Chr8 pos. 22020294 (intronic), Chr10 pos. 81316603 (3’ UTR), and Chr10 pos. 81353921 (Intergenic)) from GP7 (Table 5) represented by the surfactant genes: SFTPC, SFTPA1, and SFTPA2, respectively. A search of these variants in the family and unrelated ARDS showed that there was a higher proportion of the unrelated ARDS cohort that had all 3 variants as compared to the family control cohort. Only 1 family member (GP4) had 2 out of the 3 variants (Chr8 pos. 22020294 (intronic) and Chr10 pos. 81353921 (intergenic)), whereas at least 5 of the unrelated ARDS individuals had all 3 variants.

Emergence of clusters of orthologous genes across ARDS cases

A series of data mining experiments were performed using the DAVID database resource. We conducted 3 separate experiments with genes from groups (B), (C), and (D). Genes from group (A) were excluded as these were proprietary and we did not have access to the raw file. We clustered genes from list (B) and (C) using the search parameters applied for high stringency similarity threshold (85%). We successfully matched 308 and 956 genes from gene lists (B) and (C), respectively, based upon the availability of experimental and curated information in the DAVID database. The genes from (B) were categorized into 44 clusters (S5 File). The genes from (C) were categorized into 92 clusters (S6 File). In addition, a significant majority of group (B) (i.e. 67%) and group (C) (i.e. 85%) was comprised of genes that had previously been experimentally verified and identified as having polymorphisms. For group (D), 104 genes were utilized for our initial database search and 101 genes were successfully matched to having DAVID ids.

Using the same stringency parameter as previously, the genes were categorized into 8 clusters (S7 File). The largest cluster for group (D) was comprised of 28% of the genes (29 genes out of 101). This cluster was annotated as having functions in transmembrane composition. Additional information from these output suggested that 71% (72) of the genes from group (D) were annotated as variants based upon previous experimental information contained in the DAVID database; 67% (i.e. 68 genes) had an identified role as polymorphisms; 42% (i.e. 42 genes) were splice variants; and 45 of these genes were annotated by DAVID as having a role in post-translational modifications based upon the inclusion criteria that the protein is post translationally modified by the attachment of either a single phosphate group, or of a complex molecule, such as 5’-phospho-DNA, through a phosphate group. A follow up assessment using the Online Mendelian Inheritance in Man (OMIM) database, suggested potential roles for myosin light chain and myosin heavy chain genes in lung-related diseases. However, many of the other genes from groups B, C, and D revealed very little curated experimental information available in the OMIM database (S8 File).


Previous genomic studies in ARDS have focused on characterization of subjects unrelated to each other in an effort to identify global genomic patterns and signatures that could be informative for diagnosis or prognosis [20,22,29,30,42]. Numerous studies have suggested that there is a delicate relationship between the gene and local environment in ARDS pathogenesis. It has also been noted that there is a physiological balancing act that plays a role in this dynamic process [3,4,7,14,20,43]. The relationship is such that the following two events must occur to give rise to the disease: First, an individual must have a genetic propensity and second their microenvironment must undergo a struggle to execute an appropriate response to local insult such as trauma, sepsis, or infectious trigger [1-6,17-19]. It is the combinatorial effect from these two events that preferentially activates the disease. We hypothesized that ARDS-affected individuals have patterns of variants in their physiological repertoire that can be tracked, and then these would complement clinical diagnosis and/or clinical monitoring. The goals here were (1) to characterize the landscape of variants within protein coding but we also studied UTR regions in ARDS using an Exome sequencing approach; (2), determine the variations in signaling pathways across ARDS; and (3) use computational approaches to explore the functional consequences of ARDS. The findings from this study showed that : (1) there are unique variants not found in unaffected individuals but shared by ARDS cases; (2) Coordinated signaling pathways shared by ARDS cases are different from those found in unaffected individuals; (3) Clustering analysis demonstrated that the GP7 case exhibited less similarity in variant expression when compared with related family members in contrast to higher similarity observed with ARDS unrelated individuals; (4) Validation of variants from ARDS cases represent opportunities for contribution to gaps in coverage within the ARDS literature (i.e. literature curated genes).

The strategy undertaken in this study was to focus on a structured landscape. The structured landscape in this case was defined by the shared variant inheritance pattern represented across ARDS cases, and it consisted of two components: private variants (variants shared amongst ARDS cases) and public variants (variants found in everyone). For ARDS, we believe that the variant inheritance pattern helps define the ARDS landscape in part. Our primary motivation was to elucidate underlying inheritance patterns of “private variants” that could be hidden within a larger cohort. Private variants consist of variants shared amongst ARDS cases but not found in the cohort of related family members. This is in contrast to public variants, which would consist of variants found broadly in both ARDS cases and ARDS-unaffected individuals. The question posed was, does a focused landscape of ARDS display a distinct variant pattern from person to person or is there a shared variant pattern? An important consideration for this study, was that there was a documented presence of ARDS based on family history. This past medical history could be tracked across a generation as demonstrated by the family pedigree from Figure 2.

The variants identified support our hypothesis that there are important biological pathways that could help distinguish ARDS cases from one another. In order to place these variants in the context of the available curated literature, it was important to parse out subsets of gene lists from the larger data. We observed that there was an increase in significant pathways predicted for the GP7 case when the search window of genes was increased by the addition of groups (B) and (C) genes to the ARDS-literature curated genes. The 3,516 variants from the primary case contributed a much more robust pool of candidate genes to utilize in performing pathway predictions. The increase in sample size (i.e. # of genes included in search criteria) contributed to an increase in detection power. This was due to the signal-to-detection ratio becoming much higher, thereby leading to observed increases in significant pathways predicted. The groups (B) and (C) genes served as representative subsets of the GP7 enriched variants. By using these subsets of variants, we captured a novel pathway prediction: liver hyperplasia and hyperproliferation, which we hadn’t previously observed when searching solely within group A. Liver hyperplasia and hyperproliferation is widely implicated in respiratory complications [44-47] and has also been documented in cases of secondary ARDS due to conditions such as acute pancreatitis [48,49]. Follow up studies would be necessary in order to determine if the genes within this pathway persist broadly across other ARDS cases. In an effort to characterize variant patterns across all ARDS cases from this study, we parsed out the group (D) genes. The group (D) genes represented an important subset because it was a much more condensed list of variants shared across ARDS cases (i.e. GP7 and all unrelated ARDS JF cases). Using group (D), we intentionally wanted to isolate any ARDS-specific pathways that may have been missed while searching the other gene lists. Group (D) assessments also enabled us to place into context the contributions of the larger subsets of variants from groups (B) and (C). The 24 variants represented by 9 unique genes from (D) offer potential candidates for validation. Additional identification of 3 surfactant genes in ARDS cases was also interesting in that it was not evident in the family cohort with the exception of 1 family member (who had 2 out of the 3 gene variants). Surfactant genes have been implicated in ARDS across multiple studies and this finding is in line with their well-documented role [31,32,50,51].

Pathway assessments demonstrated that variants highly enriched in GP7, offered contributions to signaling pathways that previously weren’t found to be statistically significant for ARDS such as C-C chemokine receptor type 3 (CCR3) signaling, phosphoribosyl pyrophosphate (PRPP) biosynthesis, inhibition of matrix metalloproteases, choline degradation, glutamine biosynthesis, actin cytoskeleton signaling and epithelial junction signaling. These molecules have documented roles in lung pathogenesis and physiological signaling [24,27,52-54]. We speculate that these signaling molecules derived from GP7 could be important for the ARDS process. Of the 2,354 genes that were highly enriched in the GP7 case, we were also able to confirm some overlap with preexisting ones from the ARDS literature based upon our assessments with gene list A. The observed increases in the pool of ARDS genes from gene lists B (~2-fold) and C (~4.8- fold) is most likely the contributing factor to the identification of these novel pathways. Our data also confirm the presence of shared pathways between ARDS literature genes (group A) and GP7 enriched gene subsets, including many implicated in lung pathogenesis such as airway pathology, acute phase response signaling, and hepatic fibrosis. These findings are in line with current information on ARDS from the available curated literature genes [26,54-59].

Interestingly, the family heatmap comparisons illustrate that the second highest frequency of variants found was in the GP1 sample, representing one of the younger family members. This finding was strengthened by hierarchical clustering analysis shown in Figure 3A, which identified GP1 as having the highest similarity in variant expression pattern to GP7. This was a particularly curious finding because GP1 was the youngest family member, and had never been diagnosed with ARDS [15,60]. Recent findings indicate that younger persons can also get ARDS [61]. However, it is unclear what the consequences are for the variant inheritance pattern we observed in GP1 as this individual has not been diagnosed with ARDS, however further studies would be important so as to better determine what this discovery means. Additionally, a closer examination of the clade structure based upon clustering analysis shown in Figure 3A, displays a vastly different image from that of Figure 3B. To quantify the differences observed in Figure 3B, 7 out of the 10 ARDS controls had a high similarity in variant expression profiles. Whereas, 3 out of the 10 ARDS controls exhibited a much different variant expression profile as demonstrated in the heatmap assessments. Thus, the contrast between Figure 3A and 3B is such that there was a higher proportion of variants shared amongst the unrelated ARDS than within the family cohort. This indicates that basal expression patterns of variants exist, and that these are preferentially shared amongst ARDS cases. This presents an opportunity for us to better understand in future studies how this finding relates to clinical outcomes or disease severity for those outlier samples (i.e. samples clustered farther away from GP7 thereby denoting less similarity). However, the clinical implications of this observation remain unclear.

In a step towards developing a more generalized model, we applied clustering analysis for relationships between the groups of ARDS genes. The assessment of Clusters of Orthologous Genes (COGS) categories yielded pre-existing experimental information about the genes. The DAVID resource provided a streamlined approach for identifying the roles of previously discovered variants from ARDS-related genes, which was of benefit in our classification of the variants. Functional annotations retrieved showed that much of the genes, such as the myosin light chain genes, from groups A, B, and C had previously been identified as containing polymorphisms. These data are in line with previously published information while also highlight new possibilities for exploratory work and data sharing with the broader scientific community [62]. ARDS studies done on families with previous generations are very few [27]. As the genes identified carry specific variants, it will be important to understand the relationship of these variants across multiple ARDS cases. Based upon our findings, it is important to validate the variants whose function is currently unknown. The OMIM database, which is the primary inclusion criteria for the OMIM ontology, had a lack of enough experimental information on genes from groups B, C and D as relates to ARDS. This represents a gap in coverage within the literature for which one could contribute to experimental knowledge by validation and deposition of the variants from this study.


Taken together, the findings from this study promote the idea that there is a coordinated effort amongst signaling processes that underlie the pathogenesis of ARDS. This combinatorial signaling behavior has been well documented for lung pathogenesis and is implicated in ARDS. The study population represents an important demographic in helping understand the disease prevalence within a specific family in comparison to the unrelated ARDS controls. That said, the frequency of individuals sampled in the family cohort though not a large number, does show that there are some family members with a high degree of shared variant expression as compared to the primary case, GP7. We speculate that more work is needed to validate these additional variants and place them in the larger context of familial ARDS cases. These in-silico discoveries suggest a potential role for private variants shared by ARDS cases to be further explored in additional studies of ARDS. The potential outcomes would contribute to efforts geared towards deeper learning about the clinical consequences of these variants and broader thinking about cases for which family history suggests an underlying propensity for inheritance of ARDS.


The authors wish to thank Ms. Anna Salzberg for her help with some of the analysis steps.

Conflict of Interest

None to declare

Figure 1: Analysis Workflow. Outline of the workflow steps for data collection starting from the sample processing level, data filtering steps including pre-processing of the raw sequence data file and culminating in the final data analysis steps of assessments with the post-processed data.

Figure 2: Family pedigree. The pedigree represents the family for our primary case, GP7. Circles denote females and squares denote males. The empty circles and square are indicating those family members whose samples were not available for the study. GP7 is a female and represents our primary case. GP4, GP8, GP6, and GP3 are males. GP1 and GP2 are females. GP1, GP2, and GP6 are the offspring of GP7 andGP3 and are also the youngest family members as they were<30 years old at the study’s inception.

Figure 3A: Heatmap of Family members and the GP7 case. The heatmap was generated using a hierarchical clustering algorithm applied to each GP sample as compared to GP7, based on the presence (colored in red) or absence (colored in green) of each individual variant. The samples clustered closest to GP7 indicate a high degree of similarity of variant occurrence. GP1 was the most similar to GP7. GP1 and GP4 form a clade with GP7. GP2 was the least similar to GP7.

Figure 3B: Heatmap of unrelated individuals with ARDS and the GP7 case. The heatmap was generated using a hierarchical clustering algorithm applied to each JF sample as compared to GP7, based on the presence (colored in red) or absence (colored in green) of each individual variant. The samples clustered closest to GP7 indicate a high degree of similarity of variant occurrence. JF14 was the most similar sample to GP7 and JF10 was the least similar.

Figure 4: Relative gene frequency for IPA assessed gene lists. The variants from the study were categorized according to their respective genes and were subsequently placed in groups. They were grouped from A-D, based upon the following criteria: (1) Group A- comprised of 155 literature genes curated from ARDs-related studies; (2) Group B- comprised of 363 variants represented by 314 exonic genes affecting stop codon, genes with frameshift mutations and non-synonymous mutations; (3) Group C- comprised of 1,376 variants represented by 973 genes (which included all genes from Group B) plus genes found within the 3' UTR, 5' UTR, and non-coding RNAs and; (4) Group D- 128 variants that were represented by 104 genes shared by GP7 and all JF control cases from the study. The y-axis represents the total count of genes that belong to each group. The x-axis denotes each of the gene lists.

Figure 5: Novel pathway function identified from variants. Liver hyperplasia/hyper proliferation was a novel function predicted for the IPA assessed groups C, B, and D (as represented in the bar graph by C, B, and D respectively), but was not a significant pathway predicted in group A, the ARDS literature gene list. The total number of genes predicted was highest in group C, which represents all genes from group B, plus genes found within the 3' UTR, 5' UTR, and non-coding RNAs.

Figure 6: Genes found to be predictive of cell death and survival functions. The bar graph illustrates the number of genes for which top molecular/cellular function was cell death and survival. Assessments were done for genes from Group A (155), Group A+B (155+314), and Group A+C (155+973). Of the 95 genes from Group A, 82 genes were found to be shared amongst groups A, B and C. Forty-one additional genes were identified from Group B + Group A assessments, and 155 additional genes were identified in Group C + Group A assessments. Each pair of bars represents the frequency of genes in the search criteria and the frequency of genes identified as having a function in cell death and survival, respectively.

SNP (Variant) Function

# of Variants (SNPs) identified from GP7



Intronic or ncRNA_intronic








Upstream or downstream
















Variant function and the # of variants classified as having those related functions from the GP7 case. Functions were categorized based upon the filtered list of variants derived from those highly enriched in GP7.

Table 1: GP7 enriched variants and their function.

+Family member relationship to patient

Sample identifier

Frequency of variants found in sample













Uncle (paternal)













* N/A

* N/A

Unrelated ARDS case



Unrelated ARDS case



Unrelated ARDS case



Unrelated ARDS case



Unrelated ARDS case



Unrelated ARDS case



Unrelated ARDS case



Unrelated ARDS case



Unrelated ARDS case



Unrelated ARDS case



+The relationship of each family member to the GP7 case is shown. Unique identifiers for family members and unrelated JF samples as well as the frequency of variants detected in each sample compared to GP7 are shown in second and third column respectively. N/A denotes sample was unavailable for assay purposes.

Table 2: Variant frequency in GP family members and unrelated ARDS cases.

Analysis group ID (from Fig 4)

Number of Genes in group

Biological Description

Sample group tag





literature genes curated from ARDS- related studies


IPA/ARDS/literature genes






Exonic genes affecting stop codon and non- synonymous








Exonic genes from Group B plus genes from 3' UTR, 5' UTR,

and non-coding RNAs


Exon/frameshift/stopgain/nonsyn/3UTR/5U TR/ncRNAexon



Genes shared by GP7

and all JF cases

All ARDS cases shared genes

Group B, C, and D represented the genes derived from this study; Group A genes were those available in the ARDS literature for all biological pathway predictions using IPA.

Table 3: IPA analysis of groups of genes, unique identifiers, and biological functions.

Gene Name

Variant location and exonic information





















































*Denotes the same gene represented by 3 variants which was also found present in group B genes (Table 3). All other genes and variants listed in the table were unique to ARDS cases (i.e. group D) and not found in any other group.

Table 4: Genes and variants shared across ARDS cases.

Variants from Surfactant Genes

Gene exonic function

Gene Info






















The column labeled “CHROM” denotes chromosomal location; the column labeled “POS” denotes the coordinate for the position of the variant within the chromosomal location.

Table 5: Surfactant genes found in GP7 case.


  1. Force ADT, Ranieri VM, Rubenfeld GD, Thompson BT, Ferguson ND, et (2012) Acute respiratory distress syndrome: the Berlin Definition. JAMA 307: 2526-2533.
  2. Ware LB, Matthay MA (2000) The acute respiratory distress N Engl J Med 342: 1334-1349.
  3. Phua J, Badia JR, Adhikari NK, Friedrich JO, Fowler RA, et al. (2009) Has mortality from acute respiratory distress syndrome decreased over time?: A systematic review. Am J Respir Crit Care Med 179: 220-227.
  4. Rubenfeld GD, Caldwell E, Peabody E, Weaver J, Martin DP, et al. (2005) Incidence and outcomes of acute lung injury. N Engl J Med 353: 1685-1693.
  5. Katzenstein AL, Bloor CM, Leibow AA (1976) Diffuse alveolar damage--the role of oxygen, shock, and related factors A review. Am J Pathol 85: 209-228.
  6. Matute-Bello G, Frevert CW, Martin TR (2008) Animal models of acute lung Am J Physiol Lung Cell Mol Physiol 295: L379-399.
  7. Wellman TJ, de Prost N, Tucci M, Winkler T, Baron RM, et al. (2016) Lung Metabolic Activation as an Early Biomarker of Acute Respiratory Distress Syndrome and Local Gene Expression Anesthesiology. 125: 992-1004.
  8. Ashbaugh DG, Bigelow DB, Petty TL, Levine BE (1967) Acute respiratory distress in adults. Lancet 2: 319-323.
  9. Janz DR, Ware LB (2013) The needle in the haystack: searching for biomarkers in acute respiratory distress syndrome. Crit Care 17:192.
  10. Dowdy DW, Eid MP, Dennison CR, Mendez-Tellez PA, Herridge MS, et al. (2006) Quality of life after acute respiratory distress syndrome: a meta-analysis. Intensive Care Med 32:1115-1124.
  11. Yehya N, Thomas NJ, Wong HR (2019) Evidence of Endotypes in Pediatric Acute Hypoxemic Respiratory Failure Caused by Sepsis. Pediatr Crit Care Med 20: 110-112.
  12. Sweeney TE, Thomas NJ, Howrylak JA, Wong HR, Rogers AJ, et al. (2018) Multicohort Analysis of Whole- Blood Gene Expression Data Does Not Form a Robust Diagnostic for Acute Respiratory Distress Syndrome. Crit Care Med 46: 244-251.
  13. Vincent JL, Sakr Y, Sprung CL, Ranieri VM, Reinhart K, et al. (2006) Sepsis in European intensive care units: results of the SOAP study. Crit Care Med 34: 344-353.
  14. Tejera P, Meyer NJ, Chen F, Feng R, Zhao Y, et al. (2012) Distinct and replicable genetic risk factors for acute respiratory distress syndrome of pulmonary or extrapulmonary origin. J Med Genet 49: 671-680.
  15. Chiumello D, Marino A (2017) ARDS onset time and prognosis: is it a turtle and rabbit race? J Thorac Dis 9: 973-975.
  16. Constantin JM, Grasso S, Chanques G, Aufort S, Futier E, et al. (2010) Lung morphology predicts response to recruitment maneuver in patients with acute respiratory distress syndrome. Crit Care Med 38:1108-1117.
  17. Alberti C, Brun-Buisson C, Goodman SV, Guidici D, Granton J, et al. (2003) Influence of systemic inflammatory response syndrome and sepsis on outcome of critically ill infected Am J Respir Crit Care Med 168: 77-84.
  18. Brun-Buisson C, Minelli C, Bertolini G, Brazzi L, Pimentel J, et al. (2004) Epidemiology and outcome of acute lung injury in European intensive care Results from the ALIVE study. Intensive Care Med 30: 51-61.
  19. Calfee CS, Janz DR, Bernard GR, May AK, Kangelaris KN, et al. (2015) Distinct molecular phenotypes of direct vs indirect ARDS in single-center and multicenter studies. Chest 147:1539-1548.
  20. Blondonnet R, Constantin JM, Sapin V, Jabaudon M (2016) A Pathophysiologic Approach to Biomarkers in Acute Respiratory Distress Syndrome. Dis Markers 2016: 3501373.
  21. Puybasset L, Cluzel P, Chao N, Slutsky AS, Coriat P, et al. (1998) A computed tomography scan assessment of regional lung volume in acute lung injury. The CT Scan ARDS Study Group. Am J Respir Crit Care Med 158: 1644-1655.
  22. Walter JM, Wilson J, Ware LB (2014) Biomarkers in acute respiratory distress syndrome: from pathobiology to improving patient care. Expert Rev Respir Med 8: 573-586.
  23. Desai KH, Tan CS, Leek JT, Maier RV, Tompkins RG, et al. (2011) Dissecting inflammatory complications in critically injured patients by within-patient gene expression changes: a longitudinal clinical genomics study. PLoS Med 8: e1001093.
  24. Matthay MA, Zemans RL (2011) The acute respiratory distress syndrome: pathogenesis and Annu Rev Pathol 6:147-163.
  25. Meduri GU, Annane D, Chrousos GP, Marik PE, Sinclair SE (2009) Activation and regulation of systemic inflammation in ARDS: rationale for prolonged glucocorticoid therapy. Chest 136: 1631-1643.
  26. Nick JA, Caceres SM, Kret JE, Poch KR, Strand M, et al. (2016) Extremes of Interferon-Stimulated Gene Expression Associate with Worse Outcomes in the Acute Respiratory Distress Syndrome. PLoS One 11: e0162490.
  27. Reilly JP, Christie JD, Meyer NJ (2017) Fifty Years of Research in ARDS. Genomic Contributions and Opportunities. Am J Respir Crit Care Med 196: 1113-1121.
  28. Ren S, Chen X, Jiang L, Zhu B, Jiang Q, Xi X (2016) Deleted in malignant brain tumors 1 protein is a potential biomarker of acute respiratory distress syndrome induced by pneumonia. Biochem Biophys Res Commun 478: 1344-1349.
  29. Shortt K, Chaudhary S, Grigoryev D, Heruth DP, Venkitachalam L, et (2014) Identification of novel single nucleotide polymorphisms associated with acute respiratory distress syndrome by exome-seq. PLoS One 9: e111953.
  30. Tennessen JA, Bigham AW, O'Connor TD, Fu W, Kenny EE, et al. (2012) Evolution and functional impact of rare coding variation from deep sequencing of human Science 337: 64- 69.
  31. Floros J, Pavlovic J (2003) Genetics of acute respiratory distress syndrome: challenges, approaches, surfactant proteins as candidate genes. Semin Respir Crit Care Med 24: 161-168.
  32. Lin Z, Pearson C, Chinchilli V, Pietschmann SM, Luo J, et (2000) Polymorphisms of human SP-A, SP- B, and SP-D genes: association of SP-B Thr131Ile with ARDS. Clin Genet 58:181-191.
  33. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, et al. (2010) The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res 20:1297-1303.
  34. Wang K, Li M, Hakonarson H (2010) ANNOVAR: functional annotation of genetic variants from high- throughput sequencing data. Nucleic Acids Res 38: e164.
  35. Kramer A, Green J, Pollard J, Jr, Tugendreich S (2014) Causal analysis approaches in Ingenuity Pathway Analysis. Bioinformatics 30: 523-530.
  36. Toby IT, Widmer J, Dyer DW (2014) Divergence of protein-coding capacity and regulation in the Bacillus cereus sensu lato group. BMC Bioinformatics 15 Suppl 11:S8.
  37. R Core Team (2012) R: A language and environment for statistical computing. 3.4.0 ed. Vienna, Austria: R Foundation for Statistical Computing.
  38. Huang da W, Sherman BT, Lempicki RA (2009) Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 4: 44-57.
  39. Huang da W, Sherman BT, Lempicki RA (2009) Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res 37:1-13.
  40. Huang da W, Sherman BT, Zheng X, Yang J, Imamichi T, et (2009) Extracting biological meaning from large gene lists with DAVID. Curr Protoc Bioinformatics Chapter 13:Unit 13 1.
  41. Huang DW, Sherman BT, Tan Q, Collins JR, Alvord WG, et al. (2007) The DAVID Gene Functional Classification Tool: a novel biological module-centric algorithm to functionally analyze large gene Genome Biol 8: R183.
  42. Ware LB, Koyama T, Zhao Z, Janz DR, Wickersham N, et al. (2013) Biomarkers of lung epithelial injury and inflammation distinguish severe sepsis patients with acute respiratory distress syndrome. Crit Care 17: R253.
  43. Matthay MA, Ware LB, Zimmerman GA (2012) The acute respiratory distress syndrome. J Clin Invest 122: 2731-2740.
  44. Bryant BH, Zenali MJ, Swanson PE, Upton MP, Yeh MM, et al. (2016) Glutamine Synthetase Immunoreactivity in Peritumoral Hyperplasia in Liver: Case Report of a Metastatic Paraganglioma With Focal Nodular Hyperplasia-Like Changes and Review of an Additional 54 Liver Masses. Am J Clin Pathol 146: 254-261.
  45. Coopersmith CM, Lowell JA, Hassan A, Howard TK (2002) Hepatocellular carcinoma in a patient with focal nodular hyperplasia. HPB (Oxford) 4:135-138.
  46. Haber M, Reuben A, Burrell M, Oliverio P, Salem RR, et al. (1995) Multiple focal nodular hyperplasia of the liver associated with hemihypertrophy and vascular malformations. Gastroenterology 108:1256-1262.
  47. Van Kessel CS, de Boer E, ten Kate FJ, Brosens LA, Veldhuis WB, et al. (2013) Focal nodular hyperplasia: hepatobiliary enhancement patterns on gadoxetic-acid contrast-enhanced MRI. Abdom Imaging 38: 490-501.
  48. Zhou MT, Chen CS, Chen BC, Zhang QY, Andersson R (2010) Acute lung injury and ARDS in acute pancreatitis: mechanisms and potential intervention. World J Gastroenterol 16: 2094-2099.
  49. Opuszynska T (1973) Effect of the brand of dietary fat on rat's liver. V. Alkaline phosphatase in serum and the liver. Rocz Panstw Zakl Hig 24: 597-603.
  50. Al-Saiedy M, Gunasekara L, Green F, Pratt R, Chiu A, et al. (2018) Surfactant Dysfunction in ARDS and Bronchiolitis is Repaired with Cyclodextrins. Mil Med 183(suppl_1): 207-215.
  51. Silveyra P, Floros J (2012) Genetic variant associations of human SP-A and SP-D with acute and chronic lung injury. Front Biosci (Landmark Ed) 17: 407-429.
  52. Gao X, Qian P, Cen D, Hong W, Peng Q, et al. (2018) Synthesis of phosphatidylcholine in rats with oleic acid- induced pulmonary edema and effect of exogenous pulmonary surfactant on its De Novo synthesis. PLoS One 13: e0193719.
  53. Hove-Jensen B, Andersen KR, Kilstrup M, Martinussen J, Switzer RL, et al. (2017) Phosphoribosyl Diphosphate (PRPP): Biosynthesis, Enzymology, Utilization, and Metabolic Microbiol Mol Biol Rev 81.
  54. Juss JK, House D, Amour A, Begg M, Herre J, et (2016) Acute Respiratory Distress Syndrome Neutrophils Have a Distinct Phenotype and Are Resistant to Phosphoinositide 3-Kinase Inhibition. Am J Respir Crit Care Med 194: 961-973.
  55. Bhargava M, Becker TL, Viken KJ, Jagtap PD, Dey S, et al. (2014) Proteomic profiles in acute respiratory distress syndrome differentiates survivors from non-survivors. PLoS One 9: e109713.
  56. Kovach MA, Stringer KA, Bunting R, Wu X, San Mateo L, et al. (2015) Microarray analysis identifies IL-1 receptor type 2 as a novel candidate biomarker in patients with acute respiratory distress syndrome. Respir Res 16: 29.
  57. Evans CR, Karnovsky A, Kovach MA, Standiford TJ, Burant CF, et al. (2014) Untargeted LC-MS metabolomics of bronchoalveolar lavage fluid differentiates acute respiratory distress syndrome from health. J Proteome Res 13: 640-649.
  58. Hemnes AR, Zhao M, West J, Newman JH, Rich S, et al. (2016) Critical Genomic Networks and Vasoreactive Variants in Idiopathic Pulmonary Arterial Hypertension. Am J Respir Crit Care Med 194: 464-475.
  59. Korrodi-Gregorio L, Soto-Cerrato V, Vitorino R, Fardilha M, Perez-Tomas R (2016) From Proteomic Analysis to Potential Therapeutic Targets: Functional Profile of Two Lung Cancer Cell Lines, A549 and SW900, Widely Studied in Pre-Clinical Research. PLoS One 11: e0165973.
  60. Pediatric Acute Lung Injury Consensus Conference G (2015) Pediatric acute respiratory distress syndrome: consensus recommendations from the Pediatric Acute Lung Injury Consensus Pediatr Crit Care Med 16: 428-439.
  61. Yehya N, Keim G, Thomas NJ (2018) Subtypes of pediatric acute respiratory distress syndrome have different predictors of mortality. Intensive Care Med 44:1230-1239.
  62. Szilagyi KL, Liu C, Zhang X, Wang T, Fortman JD, et al. (2017) Epigenetic contribution of the myosin light chain kinase gene to the risk for acute respiratory distress syndrome. Transl Res 180: 12- 21.

Citation: Toby IT, Thomas NJ, Thorenoor N, Spear D, DiAngelo S, et al. (2020) Characterizing a Focused Landscape of Familial Acute Respiratory Distress Syndrome. Biomark Applic 04:141. DOI: 10.29011/2576-9588.100041

free instagram followers instagram takipçi hilesi