s2fieldsofstudy: Biology Psychology

text:
Imaging Genetic Based Mediation Analysis for Human Cognition

The brain connectome maps the structural and functional connectivity that forms an important neurobiological basis for the analysis of human cognitive traits while the genetic predisposition and our cognition ability are frequently found in close association. The issue of how genetic architecture and brain connectome causally affect human behaviors remains unknown. To seek for the potential causal relationship, in this paper, we carried out the causal pathway analysis from single nucleotide polymorphism (SNP) data to four common human cognitive traits, mediated by the brain connectome. Specifically, we selected 942 SNPs that are significantly associated with the brain connectome, and then estimated the direct and indirect effect on the human traits for each SNP. We found out that a majority of the selected SNPs have significant direct effects on human traits and discussed the trait-related brain regions and their implications.

INTRODUCTION
Advances in genetics help identify the genetic contribution to human cognition and understand the influence of specific genetic variants on cognition. The analysis of genetic variations, e.g., single nucleotide polymorphisms (SNPs), is to thoroughly understand phenotypic characteristics and the genetic mechanism of both normal and disordered brain function and behavior (Shen and Thompson, 2019). However, as the cognitive phenotypes are usually self-reported and not explicit, potentially weakening the genetic effect (Fossella et al., 2002;Goldberg and Weinberger, 2004;Green et al., 2008), the mechanisms regarding how genetic variations affect such cognitive phenotypes are not clear yet (Bi et al., 2017).
To bridge the gap between genetic variations and cognitive phenotypes, neuroimaging genetics has attracted dramatic attention (Elliott et al., 2018). In particular, with the rapid development of brain imaging acquisition techniques, a number of large-scale imaging genetic databases have been established, including the Human Connectome project (HCP) study (Van Essen et al., 2013), the UK biobank (UKbb) study (Miller et al., 2016), and the Alzheimer's Disease Neuroimaging Initiative (ADNI) study (Mueller et al., 2005), among many others. By utilizing the multi-scale data available in these imaging genetic databases, genetic variations are integrated with multimodal brain imaging, coupled with clinical and environmental factors to investigate how genes are expressed through the neuroimaging based measures and identify genetic contributors to brain activities and structures associated with cognition or neurological disorders (Nathoo et al., 2019).
Recently, some neuroimaging genetics studies have been witnessed focusing on the intersection of cognitive neuroscience and behavioral genetics (Green et al., 2008(Green et al., , 2013Bi et al., 2017;Luo et al., 2018), in which the goal is to specify the pathway that genetic effects on cognitive phenotypes are mediated by specific brain functions (Green et al., 2013). To achieve the goal, some statistical analysis approaches have been developed, typically consisting of two key components: (i) genome-wide association studies (GWAS) to identify genetic variants that have specific effects on specific brain areas; and (ii) mediation analysis to test if those brain areas produce a cognitive function (Green et al., 2013). Guen et al. (2018) explored the effect of genes on cognitive ability and found that brain activation patterns for cognitive traits can be genetic and also some regions' activation in the brain and cognitive abilities share the same genetic characteristics. He et al. (2021) proposed the possible linkage and relationship between genetic variants, brain morphometry, and working memory performance and found an SNP which might be influential on working memory. However, few works have been conducted to test the neuroimaging mediation in terms of the brain connectome although the study of brain networks becomes increasingly significant (Fornito et al., 2015;Lynn and Bassett, 2019). The main challenge comes from the construction of efficient representations for the brain connectome (Park and Friston, 2013). Typically, the functional connectivity can be measured from the statistical dependencies between physiological measures of brain activity (Van Den Heuvel and Pol, 2010) while the structural connectivity can be measured through diffusion tractography (Kazumata et al., 2019). However, the constructed connectivity graph is critically sensitive to the choice of brain atlas, including the number of region of interest (ROI) and their definitions (Balsters et al., 2016). Furthermore, unlike the voxel-wised or region-based neuroimaging phenotypes, it is difficult to conduct the GWAS for brain connectivity graph-based phenotypes to identify the corresponding potential genetic variants (Nathoo et al., 2019).
To address these issues, in this paper, we develop an imaging genetic-based mediation analysis framework to investigate the genetic effects on human cognition phenotypes mediated by an efficient representation of the brain connectome. Specifically, we focus on the HCP neuroimaging genetic study which facilitated many advancements in the field of brain networks (Craddock et al., 2015). To build up the representation of the brain connectome, an unsupervised statistical learning approach is proposed based on the tractography results extracted from the diffusion magnetic resonance imaging (MRI) data. GWAS are conducted for each brain connectome phenotype to detect the significant SNP-connectome pairs. For selected causal SNP, a regression analysis is conducted to test the mediation of the tied brain connectome via checking if it significantly affects the cognitive phenotypes.
The proposed method brings three contributions. First, a data-driven based representation of the brain connectome is developed, which possesses potential power in human cognition prediction. Second, the efficient representation can help detect genetic signals that are neglected. Third, the proposed framework can also be applied to other neurological disorders related to neuroimaging genetic databases, e.g., ADNI study, to understand the corresponding gene-connectome-cognition pathway. In conclusion, we provide a framework to study the casual pathway from SNPs data to cognitive traits, with the brain connectome as the mediator. Our result will identify significant mediation pathways, which can shed light on future investigation of SNPs→connectome→cognition mechanics. We will suggest SNP candidates that can affect human cognitive ability through altering brain structure, which will benefit further in future studies. The structure of the paper is as follows. In Section 2, we will describe the construction of the human brain connectome and the method used for GWAS and the linear modeling of effects. In Section 3, we will describe the analysis performed on the MR dataset and genetic dataset. In Section 4, the results derived from the above analysis will be discussed. Further discussions on the results and conclusion will be included in Section 5.

METHOD
In this section, we initially develop the human brain connectome based on fiber clustering on the diffusion MRI image data. After developing the brain connectome, we perform GWAS based on the structure connectivity representation to first reduce the dimension of the genetic data and then detect a significant relationship between SNPs and brain connectome. After the selection of the significant SNPs with their paired regions, we try to explore the relationship between human cognitive traits and SNPs mediated by the brain connectome.

Preprocessing
Many studies have examined parcellation in a supervised way with an existing atlas (Guevara et al., 2012;Jin et al., 2014;Gupta et al., 2017). In Liu et al. (2021), an unsupervised populationlevel approach is developed. Borrowing from this idea, we collect all the fibers together as a large fiber set after fiber tracking, denoted as F . Let F = {f ij , k = 1, . . . , m i , i = 1, . . . , n}, where f ij denotes the j-th fiber of the i-th individual, m i denotes the total number of fibers in i-th subject, and n denotes the number of subjects. Based on the length of each individual fiber, the number of points representation within a fiber varies from 4 to 60. Most clustering methods thus involve up/down-sampling of the fibers to uniform the number of points (Guevara et al., 2012). For simplicity, we take two end points {a ij } and {b ij } of each fiber f ij as the representation. More specifically, {a ij } and {b ij } are the 3D coordinates of the two end points, respectively, which are stacked into a single 6 × 1 vector (a T ij , b T ij ) T , with the order determined by the value of their coordinates. The fibers of the i-th subject {f ij , j = 1, . . . , m i } is then represented as a matrix c i with size of 6×m i . The m i represents the number of all the fibers in i th subject, which does not have to be the same across all subjects. The individual c i s are then collected by simply horizontally stacking the matrices c i s together as C = (c 1 , . . . , c n ) with the dimension of 6 × n i=1 m i , which is the reduced representation of F .

Fiber Clustering
At this end, all the fibers are represented in C. Since the original K-means algorithm has a high computational burden, we adopt a mini-batch K-means algorithm (Cho and An, 2014). In each iteration, a random set of fibers were selected to update the centroids, until we provide a stable partition on C as A K = {A (1) K , . . . , A (K) K }, where K is the number of cluster and each A (k) K interpreted as the k-th cluster. In the brain connectome studies, the parcellation number K varies from two to sometimes tens of thousands. Currently, the choice of exact cluster number remains an open problem (Eickhoff et al., 2015). A discussion about the choice of cluster number was found in Liu et al. (2021), which provides a guideline to reduce the prediction error on some specific human traits while preserving stable connectivity features.

Brain Connectome
After the fiber clusters are acquired, a subject-wise connectome is defined based on the distribution of the fibers across population clusters. Specifically, the connectome for the i-th subject is defined as a K dimension vector ω i , where K is the number of clusters. For k = 1, . . . , K, define ω ik as the proportion of number of the i-th subject's fibers within cluster k to the total number of fibers of the i-th subject, The connectome of subject i is then defined as ω i = {ω i1 , . . . , ω iK }. One of the interpretations about the connectome with K clusters is that they are distinct groups of associations pathways on the granularity level induced in our brain. With this interpretation, ω i seeks to comprehensively capture the proportion of these individual associations in i-th subject's brain.
In addition, it is important to notice that a larger ω ik can have two non-exclusive implications: (i) i-th subject has a denser axons presence in cluster k compared with other clusters; and (ii) i-th subject has a bigger volume of the relevant pair of gray matter areas that cover more fibers. This differs from our connectome method with existing functional connection based methods.

GWAS Based on Connectome Representation
First, the potential causal SNPs are identified such that the dimension of the genetic data can be reduced from a very large scale to a moderate scale and then the significant connectivity-SNP pairs are detected. Let G be the set of N G single nucleotide polymorphisms (SNPs). For the i-th subject, let x g i take the values of 0 (no minor allele), 1 (one minor allele), and 2 (two minor alleles), indicating the genetic data at the g-th locus in G, g = 1, . . . , N G , and Z i = (z i1 , . . . , z ip ) ⊤ be a p × 1 vector including the clinical confounders, e.g., age and gender. To adjust for population structures, we also included the scores for the top 2 principal components as the confounders.
Based on image-based human brain connectome, assume that K brain parcellations are derived from the diffusion MRI scans and ω i = (ω i1 , . . . , ω iK ) ⊤ is the structure connectivity representation for the i-th subject, i.e., a K × 1 vector containing the individual cluster weights satisfying ω ik ≥ 0 and k ω ik = 1. The next step is to establish a linear relationship between each SNP x g i and connectome ω ik .
where i = 1, · · · , n, k = 1, · · · , K, and g = 1, · · · , N G . Z i denotes confounders, β 0 is intercept and β k , γ k are the coefficients in the above equation. ǫ ik is the error and we assume that where ǫ ik s are iid from normal distribution with constant variance σ 2 1 . We have a huge number of genes out of which a small number of genes might have effect on ω ik hence in the model not all SNPs x g i have significant effect on ω ik . To have an initial screening, we can test the hypothesis of H 10 : β k = 0 vs. H 11 : β k = 0. By testing the above hypothesis, we can reduce the number of relevant genes in our model and hence reduce the dimension by dropping the SNPs where the p-value is large for the corresponding coefficients. Dimension reduction by dropping insignificant SNPs can help us increase the computation efficiency and it also makes the estimates biased and power gets reduced.
The subjects for which both SNP data and connectome data are available have been taken into consideration for the analysis. The estimation of the parameters in regression Equation (1) was carried out in open-source software PLINK (https://www. cog-genomics.org/plink) that is the whole genome data analysis toolset and cluster-wise results were derived containing the estimates of the parameter β k and its corresponding p-value. The selection of SNPs was made by testing the hypothesis mentioned above for parameter in Equation (1) and discarding those SNPs for which the hypothesis has been accepted. To identify significant variants, we used false discovery rate (FDR) control to select potential SNPs from our GWAS results. In this step, 942 SNPs were selected controlling for an expected FDR of 0.05, and the FDR controls were done for each cluster separately. Since the SNPs selected in this way can have a correlation with each other, we also fed our GWAS result into FUMA (Watanabe et al., 2017a) to identify independent lead SNPs with a measure of Linkage disequilibrium r 2 ≤ 0.1. The top SNPs that appear both in our FDR-selected SNPs list and in the lead SNPs of FUMA are reported in Table 1.

Genome-Wide Linear Modeling of Effects on Human Traits
In the above subsection, we have selected SNPs based on different brain regions ω ik . By matching the correspondence, we can find the related brain region with respect to each SNP. For example, for subject i = 1, suppose we find pairs of ω 11 with x 11 , x 12 ; ω 12 with x 12 , x 13 ; and ω 13 with x 14 , x 11 , then we can get a conclusion of the pairs of x 11 with ω 11 , ω 13 ; x 12 with ω 11 , ω 12 ; x 13 with ω 12 ; and x 14 with ω 13 . Thus for each significant SNP x g i , g ∈ G * and G * ⊂ G, we have its corresponding active connectome label set S g ⊂ {1, 2, · · · , K}, which will be used in the model in the next step. We adopted the additive minor allele coding to transform SNPs' information into discrete covariates. For example, if the minor allele is G, then the SNP "AA" would be encoded as 0; similarly, SNPs "AG" and "GG" would be coded as 1 and 2. We used PLINK's built in functionality to achieve this and carried out subsequent analysis. As shown in Figure 1, the brain connectome is considered a mediator in our model, which transmits the effect of genetic exposures to human traits (Mackinnon et al., 2007). To complete the rest of the model, we build a linear relationship between human traits Y, SNPs X, brain connectome ω, and confounder Z. Here, the regression is based on paired SNPs x g i and brain connectome ω ik . Brain structure-related confounding factors Z i : gender and age are also assumed to have an direct impact on human traits.
For the i-th subject, the genome connectome linear model for human trait is where i = 1, · · · , n, k = 1, · · · , K, and g ∈ G * . Z i denotes confounders and the α k are the coefficients for the active connectome ω ik . η 0 is the intercept, η and ξ are the parameters while ǫ i is the error in the above equation with the assumption that where ǫ i s are iid from a normal distribution with constant variance σ 2 2 . We estimate the parameters in the above Equation (2) by minimizing the error sum of squares. In Equation (2), we test the hypothesis H 0 : η = 0 vs. H 1 : η = 0 to check whether the parameter η is significant or not and we get the p-value of the estimate of η. For the significance of indirect pathway, we test the hypothesis H 10 : β k = 0 vs. H 11 : β k = 0 in Equation (1) and H 20 : α k = 0 vs. H 21 : α k = 0 in Equation (2). The adjusted pvalues by FDR correction are obtained with the corresponding point estimates of those parameters.

DATA ANALYSIS
In this paper, we used the subjects from the HCP 1200 Subjects Release (S1200) (https://www.humanconnectome.org/ storage/app/media/documentation/s1200/HCP_S1200_Release_ Reference_Manual.pdf) that have all four categories of data (genetic exposures, diffusion MRI, confounders (gender and age), and cognitive traits). We used the diffusion MRI to build up the population-wise brain connectome, and combine it with the rest of the data to carry out the mediation analysis. The HCP data contains 298 (149 pairs) of genetically confirmed MZ (Monozygotic) twins and 188 subjects (94 pairs) of genetically confirmed DZ (dizygotic) twins. The GWAS studies would drop one of the MZ twins (Lowe et al., 2009;Parsons et al., 2013) to prevent such twin structure to increase the type I error rate. We followed this practice and kept only one randomly chosen sample from the MZ twins in our dataset. After matching the ID with imaging data, we finally have 870 individuals (160 DZ individuals, 138 one-of-MZ individuals, and 572 non-twin individuals) as our total sample size.

MR Image Dataset
The MR and behavioral dataset used in this paper is from the HCP 1200 Subjects Release. This dataset is comprised of 1,206 behavioral and MRI data from 1,206 healthy young adult participants, collected from August 2012 to October 2015. The MR data includes structural MRI, task functional Magnetic Resonance Imaging (fMRI), resting state fMRI, and diffusion Magnetic Resonance Imaging (dMRI). We focused on the diffusion MRI data and the behavioral data. The behavioral dataset contains unit tests targeting varies domains related to human behavior, which relates to alertness, cognition, emotion, motor, personality, psychiatric, and life function. A majority of these tests are developed from the (National Institutes of Health) NIH health toolbox. It can be referred to the release manual of the S1200 dataset (https://www.humanconnectome.org/study/ hc$p$-young-adult/document/1200-subjects-data-release) and the NIH toolbox's website (https://www.healthmeasures.net/ explore-measurement-systems/nih-toolbox) for the details.

Genetic Dataset
The genetic dataset (https://www.ncbi.nlm.nih.gov/projects/gap/ cgi-bin/study.cgi?study_id=phs001364.v1.p1) from HCP was used for our analysis. In this paper, some data quality control operations were performed on the data: (1) subjects where more than 10% of the genotypes missing were removed; (2) variants where the missing genotype rate is greater than 10% were also removed from the data; and (3) variants that failed the Hardy-Weinberg test at the 10 −7 level of significance were also removed. Finally, a total of 2062590 SNPs are obtained after the quality control.
FIGURE 1 | Casual pathway diagram. Variables are colored black while parameters are colored gray. The independent variable X denotes the genetic exposure. The dependent variable Y is the cognitive trait. The brain connectomes ω act as the mediator, which transmits the indirect effect of X to Y . Both ω and Y are affected by confounders Z (e.g., gender, age). Source corresponding to letter Z, top-right: https://unsplash.com/photos/KhStXRVhfog; source corresponding to letter X, bottom-left: https://pixabay.com/illustrations/dna-biology-medicine-gene-163466/; Source corresponding to letter Y, bottom-right: https://www.pexels.com/photo/ close-up-photo-of-a-child-solving-a-rubik-s-cube-8471912/.

GWAS and Paired Brain Connectome
In our GWAS analysis, a total of 942 SNPs have been found significantly associated with 5 out of the 10 clusters (c1,c3,c4,c7,c8). In Table 1, the top SNPs with the lowest p-values for the major four clusters are shown. The clusters and brain images belong to a randomly selected subject and the clusters are colored arbitrarily to help distinguish them. As mentioned in Section 2.2, these top SNPs are independent of each other at r 2 ≤ 0.1.
In Figure 2, we displayed the 10 clusters of two different subjects, with each row corresponding to each subject. The overall regions for those clusters are similar but exist heterogeneity in different subjects. The different shapes of the white matter show that the second subject in the second row has a slightly skinnier skull. Also, the second subject has slightly less amount of fibers distributed in the second (green) cluster, compared to other clusters. It is worth noticing that in the definition of ω ik (Section 2.1.3), the denominator is the total number of fibers of the i-th subject. Thus, we are more interested in the individual's distribution of fibers, rather than the absolute number of fibers.
To help facilitate the understanding of the potential causal effect of the SNPs, we also used the Functional Mapping and Annotation (FUMA) for functional annotation (Watanabe et al., 2017b). We included SNPs in LD with our SNPs from the 1,000 Genomes reference panel provided by FUMA (https://fuma.ctglab.nl). We reported the functional annotation for the top SNPs corresponding to each cluster in Supplementary Table 1. We also used FUMA's gene-based test tool MAGMA (de Leeuw et al., 2015) to identify significant genes.

Causal Pathway Analysis
Mediating variables are used to transmit the effect of an independent variable into the response variable. Mediation analysis is a prominent method to estimate causal relationship, which is a way to understand the effect of a third variable on the relationship between two variables. Here, we focus on the effect of the human brain connectome on the relationship between human traits and SNPs. With an aim to estimate the average direct effect (ADE) of SNPs on human cognition ability, we test the hypothesis H 0 : η = 0 vs. H 1 : η = 0 for the parameter η in the Equation (2). We consider the p-values and if the p-values are less than 0.05, then we reject the null hypothesis of no direct effect at 95% level of confidence. If the null hypothesis is rejected for an SNP, then it indicates that this SNP has a direct effect on human cognition. Now, we have selected those SNPs that have a significant effect on clusters of the human brain connectome (refer to Equation (1)). Hence, the parameter corresponding to ω ik in Equation (2) together with the parameter β k from Equation (1) determines the extent of an indirect effect of the gene on human traits. One way to estimate the significance of the indirect effect is through testing hypothesis. Testing the hypothesis H 20 : α k = 0 vs. H 21 : α k = 0 and H 30 : α k β k = 0 vs. H 31 : α k β k = 0 would FIGURE 2 | Visualization of the fiber clusters of two randomly chosen subjects. Rows 1 and 3 show clusters 1 through 10 for subject S1, while rows 2 and 4 show clusters 1 through 10 for subject S2. It shows the heterogeneity of their brains and connectomes. For our model, we have considered four human traits, i.e., oral reading score (ReadEng), list sorting score (ListSort), card sort score (CardSort), and receptive vocabulary score (PicVocab). To estimate the average indirect effect (AIE) for SNPs, we consider the parameter β k estimated from Equation (1) and α estimated from Equation (2) and add the product of the coefficients ( k∈S g α k β k ) to get an overall indirect effect. For each of the four traits, we examined the effects of the SNPs selected based on the p-value of the ADE. As shown in Table 2, we also show the percentages of the intermediate effect calculated by: abs(AIE) abs(ADE)+abs(AIE) in percentage as mentioned by Mackinnon et al. (2007). Alwin and Hauser (1975) discussed that there might be situations where the coefficients of direct effect and indirect effect will have opposite signs and hence they will counteract each other. In that scenario, ratio of indirect effect to total effect might be negative or greater than one. In this scenario, the total effect is less than the total of absolute effects and a possible solution to bypass this problem is to use the absolute value of direct effect and indirect effect. A brief summary of SNPs related to direct effect and indirect effect for the above human traits is given in Table 2. In Figure 3, an example of a significant SNP and the effect of mediation variables is shown to present the fact that the human trait PicVocab is affected by the SNP through the brain connectome.

DISCUSSION
Our brain is divided into functionally specialized regions and any cognitive task is usually involved in collaboration among two or more specialized regions (Eickhoff et al., 2015). A common assumption is that the white matter tracts among these regions provide a structural basis for the communication underlying such tasks. To check the overlaps of our clusters (Ref to Figure 2) with existing anatomically meaningful fiber bundles, we surveyed a list of ROIs that have appeared in the existing literature (Friederici and Gierhan, 2013;Gupta et al., 2017;Shin et al., 2019), which include inferior longitudinal fasciculus (ILF); corpus callosum (CC); corticospinal tract (CST); uncinate fasciculus (UF); inferior fronto-occipital fasciculus (IFOF); superior longitudinal fasciculus (SLF); arcuate fasciculus (AF). Some clusters have overlaps with only one of the ROIs in the list above. For example, clusters 2 and 10 contain only the ILF, while cluster 5 contains the CST. Other clusters might have overlaps with multiple ROIs. For example, clusters 3 and 5 contain CC and UF; clusters 8 and 9 contain multiple important language related ROIs such as SLF, AF, and IFOF (Friederici and Gierhan, 2013). Although not overlapping with any ROIs in the list, cluster 4 contains the majority of the fibers that originated from the left parietal lobe and occipital lobe. The parietal lobe can be divided into 4 subregions, among which there is the inferior parietal lobule. Inferior parietal lobule has been linked with various cognitive functions (Numssen et al., 2020), some of which are believed to be lateralized to the left hemisphere, for example, our language function (Friederici, 2016). Others have also observed the left inferior parietal lobule's importance in mathematical reasoning (Eliez et al., 2001) and perspective tasks (Arora et al., 2015). The mediation analysis is of great scientific interest to explore the extent of dependence on a mediator. Here, we have performed the mediation analysis to trace the causal pathway between SNPs and human traits through the human brain connectome. The SNPs that have both significant direct effect (DE) and indirect effect (IE) were selected according to the p-values for both effects were smaller than 0.05. In total, we found 2 and 453 SNPs (out of 942 total SNPs) for human traits ReadEng and PicVocab, respectively. The top SNPs for ReadEng and 10 SNPs for PicVocab are reported in Table 2. In Table 2, the p-values of β k (see Equation (1)) and α k (see Equation (2)) are shown to represent significant pathways through the connectome. The pair (β k , α k ) signifies a pathway through the connectome and the pair (β l , α l ) signifies another pathway through the connectome. For those SNPs with significant direct effect and indirect effect, we find that a moderate percentage (8-15%) of the global effect is mediated by the brain connectome. Among the mediation paths we found in our analysis, most of the SNPs only have mediation effects through a single cluster and the mediation through cluster 7 is of most dominant. However, this may be related to the fact that cluster 7 has the most associated (606 out of the 942 selected) SNPs in our GWAS analysis. Further analysis of other imaging genetic studies with large sample size might uncover the mediation pathways through other clusters and delineate the mediation structure through multiple pathways.
To better explain the findings of our mediation analysis, we have discussed all three possible cases as outcomes of mediation analysis: significant direct effect with insignificant indirect effect, significant indirect effect with insignificant direct effect, and significant direct effect as well as indirect effect. As our focus lies on the SNPs which have significant relationship with connectome, i.e., we have already selected the SNPs where the hypothesis H 10 : β k = 0 vs. H 11 : β k = 0 is rejected, we proceed further in the following manner. In the first case, the hypothesis H 0 : η = 0 vs. H 1 : η = 0 is rejected while we fail to reject the null hypothesis H 20 : α k = 0 vs. H 21 : α k = 0 indicating that we have significant direct effect without significant indirect effect. In the second case, we fail to reject the hypothesis H 0 : η = 0 vs. H 1 : η = 0 while the null hypothesis H 20 : α k = 0 vs. H 21 : α k = 0 is rejected indicating that we have an insignificant direct effect but we find the mediation pathway through the connectome to be significant. In the third case, we are able to reject both hypotheses indicating that the direct effect is significant as well as the indirect effect, i.e., the mediation pathway through connectome is significant. In other words, the total effect of x i , i.e., SNP is transmitted to Y i through two pathways, direct pathway (direct effect is denoted by η) and indirect pathway through connectome (x i → ω ik → Y i ) whose effect is observed as α k β k termed as indirect effect. Our main focus lies on case three, i.e., both hypotheses are rejected indicating significant direct effect and a significant indirect effect. Table 3 represents the first two cases where we see only one of the DE and IDE is significant while the full results are available in Supplementary Table 2.
To understand the causal mediation of SNPs, we selected a candidate SNP rs138291276 from the third case, i.e., SNPs with both significant direct effects and indirect effects, as shown in Table 2, this SNP has significant mediation effects on both traits: PicVocab and ReadEng. On the one hand, we found that rs138291276 is mapped to gene RPH3A, which is involved in the stabilization of N-methyl-D-aspartate receptor (NMDARs) (Franchini et al., 2019). This proposed that rs138291276's observed direct effect (DE) might be related to its expression through RPH3A. On the other hand, we also observed a moderate mediation (∼14%), which suggests an indirect pathway of this SNP on human traits. The existence of both direct and indirect pathways suggests investigation of the causal mechanics of these two paths, their potential homogeneity, heterogeneity, and interaction.
For the comprehensiveness of our discussion, we also selected a candidate SNP rs4888388, from the second case, i.e., SNPs that have a significant indirect effect but not a significant direct effect. As shown in Table 1, SNP rs4888388 (imm_16_73951649) is mapped to gene CDFP1 and associated with cluster3. In Messina et al. (2017), CDFP1 is linked with microcephaly primary hereditary (MCPH), a disease that can cause a reduction in brain size and head circumference at birth. Worst cases of MCPH often result in mild to severe mental retardation (Woods et al., 2005). This evidence suggests that gene CDFP1 may be the potential medium of the mediation pathway, from SNP rs4888388 to human connectome cluster 3, and to human cognition. As mentioned before, SNP rs4888388 has an insignificant direct effect. This may propose that most of rs4888388's ability to affect human traits, is mediated through the human connectome. Lastly, the existence of this mediation pathway within the HCP1200 dataset also provides that the CDFP1's effect on brain structure and cognition not only exists in diseased patient but also among healthy young adults.
We also tried non-parametric bootstrap to obtain the mediation effects and p-values. We used the "mediation" (https://cran.r-project.org/web/packages/mediation/mediation. pdf) package in R to get the bootstrapped estimates and pvalues. We used 500 simulations with a sample size of 870 (the same as our original sample size). Then, the FDR controls are adopted to adjust p-values for addressing the multiple test issue. We see that the bootstrapped estimates are almost identical to our original estimates while the p-values support the significance of the estimates. We have also considered a 2-fold cross-validation. Specifically, we have randomly divided the data of 870 subjects into two equal groups and fitted the regression model in Equation (2) for all selected 942 SNPs.
We found that for PicVocab, no SNP has a significant direct effect on both groups; for ReadEng, we have found three SNPs have a significant direct effect in both groups but an insignificant indirect effect for both groups. Since the reduced power, we cannot find significant pathways for cross-validation. In future work, an external data set with a larger size will be considered.
To investigate the effect of genetically similar subjects from the same family, we have randomly selected one subject from each family and sample size of 434 individuals were considered for further analysis. As discussed in Section 2.2, after using FDR, we selected 2,433 SNPs for further mediation pathway analysis. We found that out of 2433 SNPs, 555 SNPs were selected in our earlier analysis of 942 SNPs where 870 subjects were considered. Among the four human traits, we only found significant mediation pathways for PicVocab. Compared to previous results of 453 mediation pathways, we also found most mediation pathways (67 mediation pathways) in PicVocab. Out of these 67 mediation pathways, we found 12 SNPs which were present in the set of 555 SNPs. Two SNPs, rs79993170 (JHU_1.246357095) and rs184384228, have also been found to be significant in our previous analysis. Another SNP rs138843033 that turns out to be significant can also be seen in our previous analysis with a p-value of direct effect as 0.06, slightly greater than our cut-off of 0.05. Thus, we conclude that SNPs rs79993170, rs184384228, and rs138843033 are considered as valid SNP candidates for mediation pathways from a conservative point of view. We have attached the results of this analysis in Supplementary Table 3. The coefficients and p-values of direct effect and the p-values of mediation pathway are included. In the table, we have also pointed out the results of the above 3 SNPs in Supplementary Table 3, which overlap with our previous analysis. We save all the results to keep the possibility of exploring mediation pathways open.
There are a couple of limitations in the findings of our current study. First, it is challenging for our proposed framework to handle the twin/family studies because it does not take into account the bias imposed by DZ twins as they share partial genes. Although removing one of the DZ twins can avoid the correlation within twins, it is inefficient due to the reduced sample size and the loss of information contained in twins. In addition, the current GWAS result is sensitive to the sampling mechanism when randomly selecting one subject from each family. In fact, we generated another four different sampling combinations of the 434 individuals, and the number of significant SNPs were 1771, 139, 2329, and 777, respectively, which varies a lot across different sampling combinations. Therefore, it is of great importance to develop advanced statistical models and address the issues in twin/family studies. Furthermore, the results need to be replicated using large samples in future studies. Second, choosing the number of clusters is still an open problem. Although the choice of 10 clusters provides us with better performance in terms of both visualization and explainability of brain connectomes, it lacks further structural justification. Thus, it will be an interesting direction to leverage brain structure and connectivity information and develop some learning techniques in finding the optimal number of clusters in brain connectomes.

AUTHOR CONTRIBUTIONS
TC, AM, and RL contributed to the conception and design of the study. TC and AM performed the analysis and wrote the first draft of the manuscript. HZ contributed to manuscript polishing. RL supervised the project, developed the method and wrote the manuscript. All authors provided critical feedback and helped shape the research, analysis, and manuscript. All authors contributed to the article and approved the submitted version.


Tumor necrosis factor receptor regulation of peripheral node addressin biosynthetic components in tumor endothelial cells

Tumor-associated tertiary lymphoid structures are ectopic lymphoid aggregates that have considerable morphological, cellular, and molecular similarity to secondary lymphoid organs, particularly lymph nodes. Tumor vessels expressing peripheral node addressin (PNAd) are hallmark features of these structures. Previous work from our laboratory demonstrated that PNAd is displayed on intratumoral vasculature of murine tumors, and its expression is controlled by the engagement of lymphotoxin-α3, secreted by effector CD8 T cells, with tumor necrosis factor receptors (TNFR) on tumor endothelial cells (TEC). The goals of the present work were: 1) to identify differences in expression of genes encoding the scaffolding proteins and glycosyl transferases associated with PNAd biosynthesis in TEC and lymph node blood endothelial cells (LN BEC); and 2) to determine which of these PNAd associated components are regulated by TNFR signaling. We found that the same genes encoding scaffolding proteins and glycosyl transferases were upregulated in PNAd+ LN BEC and PNAd+ TEC relative to their PNAdneg counterparts. The lower level of PNAd expression on TEC vs LN BEC was associated with relatively lower expression of these genes, particularly the carbohydrate sulfotransferase Chst4. Loss of PNAd on TEC in the absence of TNFR signaling was associated with lack of upregulation of these same genes. A small subset of PNAd+ TEC remaining in the absence of TNFR signaling showed normal upregulation of a subset of these genes, but reduced upregulation of genes encoding the scaffolding proteins podocalyxin and nepmucin, and carbohydrate sulfotransferase Chst2. Lastly, we found that checkpoint immunotherapy augmented both the fraction of TEC expressing PNAd and their surface level of this ligand. This work points to strong similarities in the regulation of PNAd expression on TEC by TNFR signaling and on LN BEC by lymphotoxin-β receptor signaling, and provides a platform for the development of novel strategies that manipulate PNAd expression on tumor vasculature as an element of cancer immunotherapy.

Introduction
Trafficking of leukocytes, including T-and B-cells, into lymphoid and inflamed non-lymphoid tissues involves sequential interactions of a set of homing receptors (HR) on leukocytes with cognate ligands on blood endothelial cells (BEC) (1). During an immune response, effector T cells acquire the ability to enter inflamed tissues by upregulating HR that bind to ligands that are upregulated on activated BEC. HR and HR ligands that are utilized for T-cell infiltration into tumors have been identified (2)(3)(4)(5)(6)(7). However, tumor endothelial cells (TEC) express HR ligands at low levels (8)(9)(10), and a number of studies have showed a correlation between the levels of HR ligands on the tumor vasculature and the numbers of intratumoral T cells (9,(11)(12)(13)(14)(15)(16). A better understanding of how HR ligands are regulated in TEC may provide significant opportunities to increase the number of intratumoral T cells.
Since naïve T cells do not generally enter peripheral tissue, it had been assumed that all intratumoral lymphocytes are effectors that differentiate in tumor-draining LN and home to the tumor thereafter. However, naïve T cells infiltrated tumors that had been genetically modified to secrete homotrimeric lymphotoxin-a (LTa 3 ) (17,18) or LIGHT (19). Similar results were also obtained through the intratumoral injection of either homeostatic chemokine CCL21 (20,21) or DCs genetically modified to express this molecule (22). Naïve T cells enter lymph nodes (LN) based on their expression of L-selectin and CCR7, which bind to peripheral node addressin (PNAd) and the chemokines CCL19/CCL21, respectively (23). These HR ligands are normally expressed on specialized LN blood vessels called high endothelial venules. However, they have been detected in a variety of human tumors, and their presence is associated with a positive clinical prognosis (24)(25)(26)(27)(28)(29)(30). Recently, our laboratory demonstrated that PNAd and CCL21 were co-expressed on a small fraction of blood vessels in murine melanomas and lung carcinomas growing in multiple anatomical locations (31). The primary source of CCL21 was cancer associated fibroblasts, and to a lesser extent, TEC. PNAd + CCL21 + vessels developed spontaneously, and supported the infiltration of naive CD8 T cells that differentiated into functional effectors after intratumoral activation and significantly delayed tumor outgrowth (31, 32). Thus, PNAd + CCL21 + intratumoral vessels contribute to anti-tumor immunity by generating a selfsustaining infiltration of naïve T cells into the tumor mass.
Previous work from this lab demonstrated that PNAd expression on murine melanoma TEC was not induced by either DC or LTbR signaling, and was instead controlled by effector CD8 T cells secreting LTa 3 , which signaled through TNF receptors (TNFR) expressed on TEC (31). Here, we demonstrate that TEC expressed key PNAd glycosyl transferases and scaffolding proteins normally found in LN BEC, albeit at lower levels. We also determined the mechanisms by which TNFR regulated their expression, and that checkpoint immunotherapy augments PNAd expression on the tumor vasculature. These findings provide significant insight into the basis for PNAd biosynthesis in TEC.

B C
A FIGURE 2 PNAd expression is significantly lower on tumor-associated endothelium than on lymph node endothelium. C57BL/6 mice were S.C. or I.P. injected with B16-OVA cells. Tumors were harvested 14 days later. LN were harvested from non-tumor-bearing mice. Resected tissues were prepared for IF microscopy (A, B) or flow cytometry (C) as described in Methods. (A) Representative images of unamplified and tyramide signal amplified LN and tumor sections stained with indicated markers. Scale bar = 100 µm. (B) Quantitative image summary data for unamplified and amplified tissue sections. PNAd percentages and pixel intensities were calculated on a CD31 + mask. Data represents two experiments, n = 8 sections per group. (C) Representative flow cytometry plot, percentages, and gMFI intensities of PNAd on CD45 neg Ter119 neg PDPN neg CD31 + endothelial cells in CD45 + depleted LN and I.P. tumors suspensions. PNAd gMFIs were calculated on cells gated above the fluorescence minus one (FMO) control. Data represents one experiment, n = 5 LN or tumors per group. (B, C) Results are mean standard deviation (SD) analyzed by unpaired Welch's t-test. ns: p>0.05, *p<0.05, **p<0.01, ***p<0.001, and ****p<0.0001. relative to LN BEC, although TEC from I.P. tumors express PNAd more frequently and at higher levels than TEC from S.C. tumors.

Lower expression of PNAd on tumor endothelial cells is associated with lower expression of both glycosyl transferases and scaffolding proteins
The lower level of PNAd on TEC might reflect lower expression of either glycosyl transferases that synthesize PNAd or scaffolding proteins that display it. We flow sorted PNAd + and PNAd neg cells from I.P. tumors and LN and evaluated expression of these molecules by quantitative PCR (qPCR). Consistent with other work (49,50), the glycosyl transferases Gcnt1, B3gnt3, Fut7, Chst2, and Chst4 were all expressed at significantly higher levels in PNAd + than PNAd neg LN BEC ( Figure 3). Similarly, all of these enzymes were expressed at higher levels in PNAd + TEC than in PNAd neg TEC. With the exception of Fut7, the expression levels of these glycosyl transferases were not significantly different between PNAd neg LN BEC and PNAd neg TEC ( Figure 3). However, Gcnt1, B3gnt3, Fut7, and Chst2 were expressed~2-7-fold lower in PNAd + TEC than in PNAd + LN BEC (Figure 3), consistent with the 6-fold lower gMFI of surface PNAd on PNAd + TEC ( Figure 2). Most significantly, Chst4 expression was~42-fold lower in PNAd + TEC than in PNAd + LN BEC. Chst4 is selectively expressed in LN high endothelial venules and the major enzyme responsible for the generation of 6-sulfo sialyl Lewis X on the Core 1 O-glycan structure of PNAd (39)(40)(41)(42).
Transcript levels for scaffolding proteins Cd34 and Emcn were comparable between LN BEC and TEC, regardless of whether they expressed PNAd ( Figure 4). Podxl was also expressed comparably by PNAd + and PNAd neg cells from either LN or I.P. tumors, but the transcript levels for this molecule were significantly less in TEC than in LN BEC. On the other hand, Glycam1 and Cd300lg were more highly expressed in PNAd + cells than in PNAd neg cells in both LN BEC and TEC, although as with Podxl, they were expressed at significantly lower levels in TEC than in LN BEC ( Figure 4). Thẽ 2-7-fold differences in expression of Podxl, Glycam1, and Cd300lg between PNAd + LN BEC and PNAd + TEC are again consistent with the 6-fold lower gMFI of surface PNAd on PNAd + TEC. Collectively, these results suggest that PNAd expression on TEC is achieved by upregulation of the same subsets of glycosyl transferases and scaffolding proteins as on LN BEC. However, the magnitude of these upregulations in TEC is significantly lower in every case, and the substantially deficient upregulation of Chst4 could also affect PNAd structure.
Tumor necrosis factor receptor signaling regulates the transcript levels of some scaffolding proteins and glycosyl transferases in tumor endothelial cells.
PNAd expression on LN BEC depends on signaling through LTbR (31, 51, 52) but not TNFR ( Figure 5). In contrast, PNAd is not detected by IF on the vasculature of I.P. tumors grown in TNFR1/2 -/mice, and its expression does not depend on LTbR signaling (31). By flow cytometry, only a small percentage of TEC from I.P. tumors grown in TNFR1/2 -/mice retained PNAd, although their surface expression level was lower than that of TEC from wild-type (WT) mice ( Figure 5). This confirms that PNAd expression in most TEC depends on TNFR signaling, while a small subset of TEC apparently depends on an alternative signaling pathway. Using these cells, we determined how the loss of TNFR signaling altered expression of genes involved in PNAd biosynthesis. As expected from Figure 5, the expression of genes encoding glycosyl transferases ( Figure 6) and scaffolding proteins Glycosyl transferases are expressed at significantly lower levels in PNAd + TEC than in LN PNAd + BEC. C57BL/6 mice were I.P. injected with B16-OVA cells and tumors were harvested 14 days after implantation. LN were harvested from non-tumor-bearing mice. Endothelial cells were sorted from CD45+ depleted LN and I.P. tumor suspensions, and the expression levels of indicated RNA transcripts were determined by quantitative PCR (qPCR), as described in Methods. Data from one experiment presented as 2 -DCT relative to Hprt, n=5 LN or I.P. tumors per group. Results are mean ± SD analyzed by unpaired Welch's t-test. **p<0.01. (Figure 7) in PNAd + and PNAd neg LN BEC was unaltered in TNFR1/2 -/mice relative to their WT counterparts. In PNAd neg TEC from I.P. tumors grown in TNFR1/2 -/mice, there was no reduction in the expression of genes encoding scaffolding proteins (Figure 7) or most glycosyl transferases, although the expression of Gcnt1 increased ( Figure 6). This suggests that the TNFR-induced expression of PNAd on TEC in WT mice is a consequence of upregulation of the genes identified in Figure 3 and 4. In the small subset of PNAd + TEC from I.P. tumors grown in TNFR1/2 -/mice, Chst2 expression was 7-fold lower, but no other biosynthetic enzyme, including Chst4, was significantly changed ( Figure 6). Also, Podxl and Cd300lg expression was lower by 2-and 6-fold, FIGURE 4 Scaffolding proteins are expressed at significantly lower levels in PNAd + TEC than in LN PNAd + BEC. C57BL/6 mice were I.P. injected with B16-OVA cells and tumors were harvested 14 days after implantation. LN were harvested from non-tumor-bearing mice. Endothelial cells were sorted from CD45+ depleted LN and I.P. tumor suspensions, and the expression levels of indicated RNA transcripts were determined by qPCR, as described in Methods. Data from one experiment presented as 2 -DCT relative to Hprt, n=5 LN or I.P. tumors per group. Results are mean ± SD analyzed by unpaired Welch's t-test. **p<0.01, all other paired bar comparisons are not significant (p>0.05) and have been omitted for clarity. Tumor necrosis factor receptor signaling controls PNAd expression in most TEC but not in LN BEC. C57BL/6 WT or TNFR1/2 -/mice were I.P. injected with B16-OVA cells and tumors were harvested 14 days later. LN were harvested from non-tumor-bearing mice. Resected tissues were prepared for flow cytometry as described in Methods. Representative flow cytometry plot, percentages, and gMFI of PNAd on CD45 neg Ter119 neg PDPN neg CD31 + endothelial cells in CD45 + depleted LN and I.P. tumors suspensions. PNAd gMFIs were calculated on cells gated above FMO control. Data represents one experiment, n = 5 LN or tumors per group. Results are mean ± SD analyzed by unpaired Welch's t-test. ns: p>0.05, *p<0.05, **p<0.01. respectively, and there was a trending but insignificant reduction in Glycam1 (Figure 7). These results are consistent with the possibility that an alternative signaling pathway, operating in the absence of TNFR signaling, upregulates most of the relevant genes identified in Figures 3 and 4, but is selectively deficient in its ability to upregulate Chst2, Podxl, and Cd300lg, and potentially Glycam1.
Checkpoint immunotherapy augments PNAd expression on I.P. tumor vasculature Checkpoint immunotherapy has also been associated with enhanced expression of HR ligands on the tumor vasculature, such as ICAM-1 and VCAM-1 (55, 56). To determine whether checkpoint immunotherapy altered the expression of PNAd on tumor vasculature, we treated I.P. tumor-bearing WT mice with either anti-PD-L1 monotherapy or the combination of anti-CTLA4 and anti-PD1 and analyzed PNAd expression after 14 days of outgrowth. Both treatments significantly increased the percentage of CD31 + TEC pixels that co-stained with PNAd ( Figures 8A, B). Also, the staining intensity of PNAd on CD31 + TEC pixels was significantly higher in treated tumors. Despite these increases, PNAd staining remained localized with tumor associated tertiary lymphoid structures (TA-TLS), and was not found on non-TA-TLS associated vasculature (57). Previously, we demonstrated that checkpoint immunotherapy augmented the number of T-cells in I.P. tumors (57). Since PNAd expression in I.P. tumors depends on effector CD8 T cells secreting LTa 3 (31), these results suggest that checkpoint immunotherapy augments PNAd expression on the tumor vasculature by increasing the representation of CD8 T cells secreting LTa 3 in tumors, which in turn enhances the expression of PNAd associated components.

Discussion
In this report, we identified differences in the expression of genes encoding scaffolding proteins and glycosyl transferases involved in biosynthesis of PNAd that could account for its lower expression level in TEC, and its dependence on TNFR signaling. We found that TEC and LN BEC upregulated the same subsets of molecules in association with PNAd expression, but the extent of upregulation was less in PNAd + TEC than in LN PNAd + BEC, consistent with the lower level expression of PNAd on tumor endothelium. The glycosyl sulfotransferase Chst4 was particularly poorly induced, likely contributing to a significant underrepresentation of mature 6-sulfo sialyl Lewis X glycans on both Core 1 and Core 2 branches. Previously, we demonstrated that PNAd expression on the tumor vasculature was controlled by a mechanism involving intratumoral effector CD8 T cells secreting LTa 3 , which signaled through TNFR on TEC (31). Here, we showed that loss of TNFR signaling led to PNAd neg cells that were indistinguishable from their WT counterparts, suggesting that TNFR signaling regulates the same PNAd associated genes as LTbR signaling in LN BEC. However, we also found a remnant population of PNAd + TEC with normally upregulated expression of most of these genes, but impaired upregulation of genes encoding the scaffolding proteins Podxl and Cd300lg and the glycosyl sulfotransferase Chst2, suggesting the operation of a 3 rd pathway for PNAd induction. Finally, we showed that the percentage of PNAd + TEC and their surface expression level were increased by checkpoint immunotherapy treatment. This work provides insight into the mechanisms regulating PNAd biosynthesis in TEC and provides a platform to enhance its expression to support a continual influx of naïve cells, sustaining antitumor immunity.
Our results provide new insights into the control of the biosynthetic pathway leading to expression of 6-sulfo-sialyl Lewis X on TEC. We found that B3gnt3, Gcnt1, Fut7, Chst2, and Chst4 were expressed at higher levels in PNAd + TEC relative to their PNAd neg counterparts. This is consistent with their elevated expression in PNAd + vs. PNAd neg LN BEC (49,50), and their essential roles in PNAd biosynthesis in these cells (33). However, while Chst2 and Chst4 were comparably expressed in PNAd + LN BEC, Chst4 was significantly under expressed relative to Chst2 in PNAd + TEC. Previous studies using knockout mice have established that Chst2 and Chst4 play complementary and partially redundant roles in PNAd expression in LN (39,40). Chst4 is primarily responsible for PNAd that is expressed on the luminal endothelial surface (39), while Chst2 contributes to PNAd expressed on both luminal and abluminal surfaces (40)(41)(42). Our previous observation of PNAd expression on both surfaces of tumor endothelium (31) is consistent with the very low level of Chst4 expression in TEC, and suggests that Chst2 is largely responsible for synthesizing 6-sulfo-sialyl Lewis X in TEC. The MECA-79 antibody detects 6-sulfo-sialyl Lewis X in the Core 1, but not Core 2, O-glycan biantennary branch (37). Since Chst2 mediates GlcNAc sulfation on the Core 1 branch inefficiently (41), the levels of 6-sulfo-sialyl Lewis X on TEC may be higher than is suggested by MECA-79. PNAd structures on TEC may also be largely "single-armed", with only Core 2 derivatized by 6-sulfosialyl Lewis X. Regardless, PNAd expression on Peyer's Patch endothelium is entirely dependent on Chst2 (39,40). This, together with our earlier work (31, 32), demonstrates that these structures support naïve lymphocyte entry.
We found that PNAd + TEC and LN BEC also express similar scaffolding proteins that have been shown to display PNAd and regulate them similarly. Cd34 and Emcn were expressed comparably by PNAd + and PNAd neg endothelial cells from both LN and I.P. tumors. In contrast, Podxl expression also did not differ between PNAd + and PNAd neg endothelial cells in either tissue, but its overall level of expression was substantially higher in LN BEC than TEC. Finally, Glycam1 and Cd300lg were expressed at higher levels in PNAd + than in PNAd neg cells in both tissues, and again, the transcript levels for these molecules were significantly less in PNAd + TEC than in PNAd + LN BEC. These scaffolding proteins are redundant with one another in PNAd display, and at least partially redundant in promoting Lselectin engagement with PNAd, as mice deficient in Cd34 (58) or Glycam1 (33) show no impaired trafficking of naïve T cells to peripheral LN. However, together with the lower expression levels of multiple glycosyl transferases, particularly Chst4, the overall lower level of Podxl, Glycam1, and Cd300lg in PNAd + TEC may also contribute to the lower level of PNAd expression on the tumor vasculature.
This report also identified components of the PNAd biosynthesis pathway in TEC that are regulated by TNFR signaling. The number of PNAd + TEC is reduced by about 80% in TNFR1/2 -/mice. The resulting PNAd neg cells showed no reduction in any evaluated glycosyl transferases or scaffolding proteins and are thus indistinguishable from their counterparts in tumors from WT mice. By comparison of gene expression in PNAd + and PNAd neg LN BEC and TEC, this suggests that TNFR signaling induces PNAd by acting on the same genes as LTbR signaling does in LN BEC, albeit less efficiently, and with a particular deficiency in upregulation of Chst4. However, a small percentage of TEC from these tumors continue to express a lower level of PNAd. These cells selectively express lower levels of Chst2, Podxl, and Cd300lg, and potentially Glycam1, all of which are also elevated in PNAd + cells vs PNAd neg cells from WT mouse tumors. It is possible that TEC that retain PNAd in the absence of TNFR signaling are responsive to signals from LTbR as an alternative. However, in PNAd + LN BEC, LTbR signaling upregulates Chst2 and Glycam1, but also Gcnt1, Chst4, Fut7, and not Podxl or Cd300lg (51, 52). This difference in response pattern seems inconsistent with the hypothesis that the remnant fraction of PNAd + TEC in TNFR1/2 -/mice are regulated by LTbR signaling. Nonetheless, these results demonstrate a distinct regulation of PNAd biosynthetic components by TNFR1/2 and LTbR signaling. the mutually exclusive responses of LN BEC and TEC to these two different signals also point to important anatomic microenvironmental controls on endothelial cell signaling responsiveness.
Related to this, the reasons that all evaluated PNAd glycosyl transferases and several scaffolding proteins are expressed at significantly lower levels in PNAd + TEC than in LN PNAd + BEC are unclear. It is possible that this reflects the relative efficiency with which the TNFR and LTbR signaling pathways engage transcription factors that upregulate these genes. On the other hand, several studies have also shown that the epigenetic profile of BEC is different from TEC (reviewed in (59). Two studies demonstrated that human umbilical vein endothelial cells cultured with conditioned tumor media underwent epigenetic modifications that resulted in reduced VCAM-1 and ICAM-1 expression (60, 61). Thus, epigenetic modifications that occur in the tumor microenvironment could also reduce upregulation of PNAd associated components in TEC. Finally, it is possible that this is controlled by the immunosuppressive microenvironment of the tumor. In murine methylcholanthrene-induced fibrosarcoma, regulatory T cells have been shown to limit PNAd expression on tumor vessels (62,63), although the mechanism of action is unknown. Conversely, checkpoint immunotherapy has been associated with enhanced expression of HR ligands, such as ICAM-1 and VCAM-1, on tumor vasculature (55, 56), and we found that anti-PD-L1 monotherapy or the combination of anti-CTLA4 and anti-PD1 enhanced the fraction of CD31 + TEC expressing PNAd and their surface levels of PNAd. However, the fraction of PNAd + CD31 + TEC and their PNAd expression level in treated tumors still did not reach that of LN BEC. In related work, we demonstrated that these therapies also led to an increase in the number and size of TA-TLS, but we found no evidence for PNAd expression on non-TA-TLS associated tumor vasculature (57). Since PNAd expression in I.P. tumors depends on effector CD8 T cells secreting LTa 3 (31), and these therapies also enhanced the representation of T-cells in I.P. tumors (57), this suggests that immunotherapy-augmented PNAd expression on the tumor vasculature is driven by an increased number of intratumoral CD8 T cell effectors secreting LTa 3 . It also suggests that the increased fraction of TEC expressing PNAd create nucleation sites for the formation of new TA-TLS, or the expansion of existing ones. Determining the factors that limit the expression of PNAd scaffolding proteins and glycosyl transferases could identify targets for enhancing the expression of these molecules to enhance naïve and central memory cell infiltration into tumors, promote the development and/or expansion of tertiary lymphoid structures, and augment anti-tumor immunity.

Materials and methods
Mice Female C57BL/6 mice were from the National Cancer Institute. TNFR1/2 -/mice were from the Jackson Laboratory. All mice were bred and maintained in specific pathogen-free conditions. All experiments were carried out on female mice that were~8-12 weeks of age. All protocols and experiments were approved by the University of Virginia Institutional Animal Care and Use Committee.

Digestion of resected tissues
Resected LN and tumors were minced and digested with a solution of 0.1 mg/ml DNase I (Sigma), 0.8 mg/ml Collagenase Dispase (Sigma), and 0.2 mg/ml Collagenase P (Sigma) for 30 minutes at 37°C. Every 5 minutes, tissue suspensions were pipetted up-and-down several times. Digested tissues were depleted of red blood cells using RBC Lysing Buffer Hybri-Max (Sigma) according to manufacturer's instructions.

Enrichment of cells from digested tissues
Digested tumor suspensions were depleted of hematopoietic cells using CD45 magnetic beads (Miltenyi Biotec) on an AutoMACS Pro Separator (Miltenyi Biotec) according to manufacturer's instructions. CD45 + depleted suspensions were stained with biotinylated anti-CD31 (0.5 µg/mL) for 15 minutes at 4°C and enriched using anti-biotin magnetic beads (Miltenyi Biotec) on an AutoMACS Pro Separator according to manufacturer's instructions. BEC and TEC were sorted to highest purity according to the procedure below.

Flow cytometry and cell sorting
Cell surface staining was done in PBS containing 2% FBS, 2 mM EDTA (Sigma), and 2 mM NaN 3 (Sigma) for 30 minutes at 4°C. Live/Dead Aqua (Invitrogen) or 4,6-diamidino-2phenylindole (Sigma) were used to exclude dead cells from analysis. Endothelial cells were defined as live, singlet, Ter119 neg , CD45 neg , PDPN neg , CD31 + . Samples were run on a FACSCanto II (BD) or Attune NxT (ThermoFisher/Invitrogen) and analyzed using FlowJo Software (BD Bioscience). For qPCR experiments, pre-enriched endothelial cell populations were sorted on an Influx Cell Sorter (BD) directly into RNAlater Stabilization Solution (ThermoFischer Scientific) or PBS. A small aliquot of the sorted population was re-run to determinẽ 95% purity.

Four-color IF microscopy
Preparation of tumor tissue for IF microscopy using formalin fixation and cutting frozen blocks on a -20°C cryostat has previously been described (65). Tumor sections were incubated in 100% methanol for 10 minutes at -20°C. Then, slides were immersed in PBS for 10 minutes at room temperature. For Fc blocking, tumor sections were incubated with 0.5 mg/mL anti-CD16/32 (BioXcell; Clone:2.4G2) unconjugated antibody in PBS containing 5% BSA (Sigma) and 0.3% Triton X-100 (Sigma) for 15 minutes at room

Four-color IF image acquisition
One section per tumor was evaluated and all images were captured on an AxioImager with Apotome (Zeiss). Fluorescence minus one (FMO) or isotype staining controls were used to establish thresholds and exposure times to visualize positive signals while minimizing background fluorescence. Fifteen-20 low magnification images per tumor section were captured and stitched together using ImageJ Software (NIH) to create an image of the entire tumor section. Values across respective regional images were averaged together for each tumor section. In all experiments, the quantification of CD31 + and PNAd + areas were performed using ImageJ software (NIH) on original fluorescence images taken at identical exposures. Consistent thresholds were applied to each image to identify CD31 + and PNAd + pixels. The percentages and gMFI of PNAd + pixels within the region of interest of CD31 + area was calculated using ImageJ Software (NIH). For image presentation, brightness and contrast were linearly adjusted and colormerged images were generated using Photoshop CS6 Software (Adobe).

Quantitative RT-PCR
RNA was purified from flow sorted endothelial cells using RNEasy kits (Qiagen). High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems) and purified RNA was used to generate cDNA. Amplification was performed using TaqMan Fast Advanced Master Mix (Applied Biosystems) and QuantStudio 6 Flex Real-Time PCR system (Applied Biosystems) with the following program: 50°C for 2 minutes; 95°C for 2 minutes; 40 cycles of 95°C for 1 second, 60°C for 20 seconds.

Statistical analyses
Statistical details of each experiment in this work are reported in the figure legends. Normality of data distribution was determined by D'Agostino-Pearson omnibus normality test and variance between groups was assessed by the F-test. P-values for the comparison between two independent groups were calculated by Welch's t-test. Error bars shown in graphical data represents mean ± standard deviation (S.D.) for normally distributed data. P<0.05 was considered statistically significant. All statistics were calculated using Graph Pad Prism version 7.0, R version 4, and the SAS software suite version 9.4.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

Ethics statement
The animal study was reviewed and approved by University of Virginia Institutional Animal Care and Use Committee.

Author contributions
Experiments involving tumor implantation, flow cytometry, endothelial cell sorting, and quantitative PCR were performed and analyzed by AR and GP. Immunofluorescence image capture and analysis were performed by AR. VE supervised all experiments associated with this project and participated in experimental design and interpretation of results. All authors contributed to the article and approved the submitted version.